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Executive  summary 

We  employed  first  principles  calculations  with  supercell  approach  to  study 
various  defects  in  AIN,  GaN,  InN,  related  wide  gap  semiconductors,  and  their  alloys.  For 
the  exchange  and  correlation  functions,  we  used  both  the  standard  local  density 
approximations  (or  generalized  gradient  approximations)  and  time  consuming  HSE 
hybrid  functional  which  can  provide  electronic  band  gaps  that  are  much  closer  to  the 
experimental  values.  The  later  is  being  used  when  more  accurate  determination  of  defect 
levels  is  needed.  In  this  program,  we  also  employed  the  FEFF  codes,  which  is  the 
software  specially  designed  for  the  simulation  of  the  x-ray  absorption  spectrum,  to  allow 
a  direct  comparison  between  our  crystal  structure  models  and  actual  measurements.  We 
have  accomplished  several  important  results,  a)  We  reconfirm  our  earlier  calculations 
that  Mg  and  Be  are  deep  acceptors  in  AIN  by  using  HSE  hybrid  functional  which  does 
not  suffer  from  band  gap  underestimation  problems,  b)  We  calculated  A1  Frenkel  pairs 
(interstitial  A1  and  A1  vacancy  pairs)  in  AIN.  Both  the  stability  and  the  structure  are  very 
similar  to  the  case  of  Ga  Frenkel  pairs  in  GaN  that  we  have  studied  last  year.  We  also 
found  that  the  barrier  for  the  pair  to  recombine  is  ~0.2  eV  which  is  very  small  (similar  to 
the  case  of  GaN),  making  the  pairs  unlikely  to  be  stable  at  room  temperature,  c)  We 
studied  the  x-ray  absorption  spectrum  of  InN  and  ln203  alloys.  Special  attentions  have 
been  paid  to  understand  the  source  of  the  features  in  the  calculated  XANES.  We  found 
that  the  XANES  spectrum  of  In  in  the  ln203  structure  is  very  different  from  that  of  In  in 
wurtzite.  We  have  also  studied  the  effects  of  interchanging  the  anions  (O  by  N  or  vice 
versa).  To  understand  the  source  of  the  changes  in  the  XANES  spectra  in  each  case,  the 
decomposed  electron  density  of  states  were  studied,  d)  We  have  studied  hydrogen  in 
ln203.  In  most  semiconductors,  H  is  amphoteric,  i.e.  can  act  as  donor  or  acceptor  to 
counteract  leading  carriers.  However,  for  ln203,  we  found  that  H,  is  exclusively  a 
shallow  donor  (similar  to  the  case  of  ZnO).  In  addition,  H  atoms  can  also  occupy 
substitutional  oxygen  sites  (Ho)  and  act  exclusively  as  donors  in  this  configuration.  Our 
findings  confirm  that  hydrogen  acts  as  a  source  of  n-type  conductivity  in  ln203.  e)  The 
elastic  constants  and  sound  velocities  in  wurtzite  GaN  and  InN  as  well  as  other  wurtzite 
semiconductors  (SiC,  ZnO,  and  CdSe)  as  functions  of  pressure  are  studied.  We  showed 
that  the  elastic  constants  are  lower  in  materials  with  higher  ionicity.  Generally,  the  sound 
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velocities  follow  the  same  trend  as  the  elastic  constants,  f)  When  ZnO  is  doped  with 
large  group  V  anions,  such  as  P,  As  or  Sb,  phase  precipitations  might  occur.  This  is 
mainly  because  of  the  large  site  mismatched  and  the  high  concentration  of  the  dopant. 
We  studied  the  electric  properties  at  the  Z113P2  and  ZnO  interface  and  found  that  a  layer 
of  hole  could  be  formed.  This  layer  could  give  the  p- type  carrier  in  the  measurements, 
g)  Recently,  experimentalists  measured  the  local  vibration  modes  of  O  doped  CdTe  and 
found  two  sets  of  frequencies,  i.e.  one  line  at  350  cm’1  and  another  set  of  two  lines  at  ~ 
1100  cm’1.  The  350  cm'1  has  been  assigned  to  O  substituted  for  Te  (Oxe)  and  the  high 
frequency  lines  (-'1100  cm-1)  have  been  assigned  to  the  vibration  associated  with  Oxe- 
V cd  complexes.  Our  calculations  showed  that  the  assignment  of  the  350  cm’1  mode  is 
correct.  However,  the  Oxe-Fcd  complexes  can  give  the  highest  frequency  of  only  -460 
cm’1,  which  is  far  smaller  than  the  observed  1100  cm’1.  Therefore,  the  1100  cm’1  mode 
can  not  come  from  Ore-  Vai  complexes,  h)  For  first  principles  computations  of  alloy 
systems,  generally  a  supercell  approach  is  used.  We  investigated  the  size  of  unit  cell  that 
can  sufficiently  give  the  reliable  results.  In  this  work,  the  quaternary  alloys  AlxGayIni.x-yP 
are  used  for  the  test.  We  found  that  at  least  the  supercell  size  with  five  cation  atoms 
(including  all  cell  shapes)  has  to  be  used  to  obtain  reasonable  results. 
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(Following  work  schedule  and  milestones  are  the  same  as  previously  submitted) 

Work  schedule: 


Month  0: 

o  Start  calculations  of  Fn  and  N,  and  /7-dopants  (especially,  Mg  and  Be)  in 
AIN  and  AlGaN  using  first  principles  supercell  approach.  Mg  and  Be  are 
potential  /7-dopant  for  AIN  and  AlGaN  alloys. 

o  Start  calculations  of  defect  complexes  between  native  defects  and  dopants 
in  GaN,  AIN,  and  related  wide-gap  materials. 

o  Investigating  various  possible  defect-complex  configurations  to  identify 
the  most  probable  ones. 

Month  6: 

o  Study  stable  defects  and  defect  complexes  in  deeper  details.  Certain 
signatures  of  them  will  be  simulated  for  comparison  with  experiments. 

o  Start  calculations  of  Frenkel  pair  formation  between  Fn  and  N,  by  directly 
calculations  of  the  energy  as  a  function  of  pair  distance  using  nudge 
elastic  band  method. 


-4- 


Milestones: 


6  month: 

o  Obtained  the  results  of  at  least  one  native  defect  in  AlGaN  alloy. 

o  Obtained  the  results  of  at  least  one  defect  complex  in  GaN  or  AIN  or  related  wide 
gap  materials. 

o  Obtained  preliminary  results  on  the  complex  between  some  native  defects  and 
/^-dopants. 


12  month: 

o  Complete  results  on  Fn  and  N,-  in  AIN  and/or  AlGaN  alloy.  Some  preliminary 
results  on  the  Fn  and  N,-  Frenkel  pair  in  AIN  will  be  obtained. 

o  Results  on  the  binding  energies  and  stabilities  of  F\-N,  Frenkel  pairs  as  a  function 
of  the  pair-distance. 

o  Signatures  for  comparison  with  experimental  measurements  of  some  interesting 
defects  and  dopants  in  AIN,  AlGaN  alloys  and  related  wide-gap  materials. 
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Research  Activities 


Computational  methods  and  codes 

To  study  the  electronic  as  well  as  other  properties  of  these  semiconductor  systems,  we 
employed  first-principles  calculation  approach.  Most  of  our  calculations  are  based  on 
density-functional  theory  and  ab  initio  projector- augmented  wave  potentials  as 
implemented  in  the  VASP  code  [1].  Generally,  either  the  local  density  approximation 
(LDA)  or  the  generalized  gradient  approximation  (GGA)  of  the  Perdew,  Burke,  and 
Emzerhof  (PBE)  [2]  are  used  for  the  exchange-correlation  functions.  These  functional 
are  known  to  greatly  underestimate  the  band  gap.  In  some  cases  where  the  more  correct 
electronic  band  gap  is  desired,  we  use  the  hybrid  functional  as  proposed  by  Heyd, 
Scuseria,  and  Emzerhof  (HSE)  [3]  which  is  much  more  computational  demanding.  A 
screening  parameter  of  ju  =  0.2  A'1  is  used  for  the  screened  non-local  exchange  [4],  The 
hybrid  functional  is  a  mixture  of  25%  Hartree-Fock  exchange  functional  and  75%  PBE 
exchange  functional.  For  the  x-ray  absorption  spectroscopy  (XAS)  simulation,  we  use 
the  ab  initio  multiple  scattering  code  called  FEFF  (version  8.2)  [5].  Because  the  FEFF 
code  cannot  optimize  the  crystal  structures,  for  the  XAS  simulation  of  local  structures 
such  as  impurities  and  alloys,  it  is  important  that  atomic  relaxations  have  to  be  performed 
using  the  first  principles  calculations  (VASP)  before  importing  the  coordinates  into  FEFF. 
Depending  on  the  nature  of  each  system  studied,  the  specific  computational  details  are 
varied  but  can  be  found  in  each  attached  publication. 
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Computational  approach  to  study  defects  in  semiconductor  crystal 

Bulk  semiconductor  crystal  is  an  ideal  case  for  performing  calculations  because  the 
system  can  be  defined  by  the  smallest  possible  unit  cell  that  repeated  periodically  in  all 
directions  by  the  three  translation  vectors.  The  highly  ordered  repetition  of  the  unit  cell 
is  an  important  factor  that  defines  the  electronic  properties  of  semiconductor  crystals. 
For  bulk  systems,  the  calculations  are  done  at  the  infinite  repetitions  of  the  unit  cell.  This 
is  highly  realistic  because  in  actual  crystal  systems  the  number  of  atoms  is  in  the  order  of 
10  .  Although,  there  is  virtually  unlimited  number  of  atoms,  in  practice  we  do  not  need 
to  repeat  the  calculation  for  every  unit  cells.  We  use  the  Bon  von  Karman  boundary 
conditions  and  perform  the  calculation  on  only  a  single  unit  cell  which  most  of  the  time 
composed  of  4  atoms  or  less.  For  defect  or  impurity  study,  the  situation  is  different. 
Majority  of  the  crystal  is  still  highly  periodic  with  a  small  number  of  defects  (or 
impurities)  scattered  in  the  crystal  interrupting  the  periodicity.  All  atoms  in  the  system 
are  no  longer  equivalent  and,  in  principle,  each  and  every  atom  has  to  be  calculated 
separately.  This  is  prohibitively  demanding  in  tenns  of  computation.  To  reduce  the 
problem  into  a  manageable  one,  the  supercell  approach  is  used.  In  the  supercell 
approach,  the  unit  cell  is  repeated  a  few  times  in  each  direction  to  create  a  larger  cell 
(called  supercell)  which  serves  as  a  new  unit  cell.  For  example,  the  72  atom  supercell  is 
a  result  of  the  repetition  of  the  4  atom  wurtzite  unit  cell  three  times  each  in  ai  and  a2 
directions  and  two  times  in  a3  direction;  leading  to  the  multiplication  factor  of  3  x  3  x  2  = 
18.  This  new  unit  cell  has  a  new  three  translation  vectors,  i.e.  3ai,  3a2  and  2a2.  To  study 
a  defect,  a  defect  is  created  in  this  supercell.  Of  course,  the  defect  is  also  repeated  as  the 
supercell  repeats;  resulting  in  an  unusual  high  defect  concentration  compared  to  reality. 
Therefore,  it  is  desirable  to  use  the  supercell  size  as  large  as  computationally  possible  to 
study  defects  that  are  generally  of  low  concentrations.  In  practice,  the  supercell  size  as 
small  as  72  or  even  36  atom  can  already  provide  reasonable  results.  In  this  study,  we 
generally  use  a  supercell  size  of  96  atoms  for  the  wurtzite  crystal  and  64  atoms  for 
zincblende  crystal.  We  used  a  set  of  non-T  k-points  sampling  and  plane  wave  cut-off  of 
400  eV  for  the  plane  wave  basis  set.  In  supercell  calculations,  the  formation  energy  of  a 
defect  X  in  charge  state  q  is  defined  as 
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E'[X^  =  Emsc^)-Etc,x(buIk)  +  naMa+qEF,  (1) 

where  Etot.sdX9)  is  the  total  energy  of  a  supercell  containing  the  defect,  Etot,sc(bulk)  is  the 
total  energy  of  the  same  supercell  but  without  any  defect.  na  is  the  number  of  the  atoms 
of  specie  a  being  removed  from  (positive  value)  or  added  to  (negative  value)  a  bulk 
supercell  to  form  defect  X.  /ua  is  the  reservoir  energy  of  specie  a  (the  chemical 
potential).  EF  is  the  Fermi  level,  referenced  to  the  valence-band  maximum  of  bulk,  Ev. 
The  chemical  potentials  in  Eq.  (1)  depend  on  the  growth  conditions  and  constrained  by 
phase  precipitation  limits. 

As  an  example  we  can  look  at  the  case  of  Mg  substitution  A1  in  AIN.  Under  Al- 
rich  condition,  the  chemical  potential  of  A1  is  set  at  the  highest  value  possible  which  is 
//ai  =  //Ai[metai]  =  formation  energy  of  solid  metal  A1  (if  //Ai  is  higher  than  this  value,  the  A1 
metal  will  be  form  on  the  sample  instead  of  AIN).  In  order  for  AIN  to  be  stable  during 
the  growth,  it  is  required  that  jua i  +  // n  =  //a in  =  heat  of  formation  of  AIN.  This  fixes  //n  = 
Aain  -  //ai-  The  chemical  potential  of  Mg  is  set  at  the  highest  value  possible,  i.e.  the 
fonnation  energy  of  metal  Mg  (//Mg  =  //Mg[metai])-  The  formation  energy  of  MgAi  in  charge 
state  q  is  defined  as 

Ef  (Mg®, )  =  EtotSC  (Mg’, )  -  EtotSC  (bulk)  +  //A1  -  //Mg  +qEF 


...(2) 


Magnesium  and  beryllium  doped  AIN 

Group-II  metals  such  as  Mg  and  Be  are  shallow  acceptors  when  substituting  for  Ga  in 
GaN.  Magnesium  has  been  used  successfully  to  dope  GaN  into  /;-typc;  leading  to 
commercial  blue  lasers  in  today  devices  [6-8],  We  have  previously  suggested  based  on 
first  principles  calculations  that  beryllium  is  a  potentially  good  p- type  dopant  in  GaN  [9], 
However,  beryllium  can  potentially  form  interstitial  (Be,,,,)  and  bind  with  the 
substitutional  acceptor  (BeGa)  turning  the  acceptor  into  a  single  donor.  Previously,  we 
suggested  that  if  one  can  overcome  interstitial  Be  problem,  beryllium  is  a  potentially 
better  acceptor  in  GaN  than  Mg  [9].  Currently  the  interests  have  shifted  to  AIN.  With  a 
larger  band  gap  of  6.2  eV,  AIN  is  a  potential  material  for  UV  optoelectronic  applications. 
However,  AIN  is  naturally  n-type  and  /;-type  AIN  is  currently  unavailable.  Similar  to 
other  wide  gap  semiconductors,  AIN  also  suffers  from  doping  asymmetry,  i.e.  easy  to  be 
doped  into  one  conducting  type  (n-type)  than  another  (p-type)  [10,  11].  There  are  also 
some  preliminary  theoretical  investigations  of  doping  in  AIN  [12-14].  Because  the  lattice 
constant  of  AIN  is  very  similar  to  that  of  GaN,  both  Mg  and  Be  are  potential  candidates 
as  for  /7-type  doping  in  AIN.  Recent  experimental  results  showed  that  Mg  is  a  deep 
acceptor  in  AIN  (with  ionization  energy  of  around  0.4  -  0.6  eV)  [15].  In  our  previous 
program,  both  Mg  and  Be  as  dopants  in  AIN  have  been  studied  by  using  DFT-LDA.  Our 
results  also  showed  that  both  Mg  and  Be  are  deep  acceptor  in  AIN  with  the  ionization 
energy  of  ~0.6  eV.  In  addition,  we  also  found  that  interstitial  Mg  and  Be  (Mg,,,,  and 
Be,,„)  have  rather  low  formation  energy  and  can  cause  spontaneous  compensation 
especially  for  the  case  of  Be.  It  is  known  that  DFT-LDA  or  DFT-GGA  largely 
underestimate  semiconductor  band  gaps.  This  could  potentially  cause  wrong 
determination  of  dopant  ionization  energies.  To  ensure  the  validity  of  the  DFT-LDA 
results,  in  this  program  we  studied  Mg  and  Be  doped  AIN  by  using  HSE  hybrid 
functional  which  gives  band  gap  in  good  agreement  with  the  experiment. 

Table  I  shows  the  lattice  parameters  and  band  gap  values  obtained  from  DFT- 
LDA,  DFT-GGA  and  HSE  in  comparison  with  experiments.  All  calculations  gave 
reasonable  lattice  parameters  with  the  experimental  values.  However,  only  the  HSE 
calculations  gave  the  electronic  band  gap  in  agreement  with  the  experimental  values  (to 
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within  10%).  As  a  result,  the  formation  energies  as  well  as  the  ionization  energy  of 

impurities  and  defects  in  AIN  calculated  by  HSE  are  expected  to  be  more  reliable. 

Table  I  The  equilibrium  lattice  parameters  of  wurtzite-AIN  ( a  and  da)  and  the  band 
gap  (Eg)  given  by  LDA  and  GGA-PBE  and  HSE  for  the  exchange- correlation 
functional.  The  experimental  values  are  taken  from  Ref.  [16]. 


LDA 

GGA-PBE 

HSE 

Exp. 

a  (A) 

3.091 

3.131 

3.102 

3.112 

c/a 

1.600 

1.605 

1.603 

1.609 

Eg  (eV) 

4.37 

4.02 

5.64 

6.20 

Figure  1  shows  the  local  atomic  structures  of  MgAi  and  BeAi  (for  the  -1  charge 
state).  All  four  Mg-N  bonds  have  almost  the  same  length  of  -2.01  A  which  is  about  5% 
larger  than  the  Al-N  bond  length  of  1.89  A.  On  the  other  hand,  all  four  Be-N  bonds 
(-1.83  A)  are  slightly  smaller  than  the  Al-N  bond  length.  This  can  be  attributed  to  the 
fact  that  Mg  atom  has  larger  radii  than  that  of  A1  atom,  making  the  surrounding  N  atoms 
relax  outward  while  Be  atom  has  smaller  radii  than  that  of  A1  atom,  making  surrounding 
N  atoms  relax  inward. 


Fig.  1  Local  atomic  geometry  of  MgAi  and  BeAi  in  wurtzite  AIN.  The  large  blue 
spheres  are  A1  atoms  and  the  small  light  blue  spheres  are  N  atoms.  The  isosurface 
charge  density  of  the  Mg  and  Be  related  state  is  also  illustrated. 


The  formation  energies  of  MgAi  and  BeAi  are  shown  in  Fig  2.  The  formation 
energies  obtained  by  HSE  calculations  are  turned  out  to  be  in  agreement  with  the  values 
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obtained  earlier  by  using  DFT-LDA  to  within  0.5  eV  (not  shown).  This  indicated  that 
our  earlier  results  were  sufficiently  reliable.  The  acceptor  ionization  energy  is  defined  as 
the  value  of  the  Fenni  level  where  the  MgAi  (or  BeAi)  in  neutral  and  -1  charge  states  have 
the  same  formation  energy.  From  HSE  calculations,  we  obtained  the  ionization  energy  of 
0.43  and  0.72  eV  for  MgAi  and  BeAi,  respectively.  These  reconfirmed  that  both  MgAi  and 
BeAi  are  deep  acceptors  in  AIN.  Our  ionization  energies  are  in  a  reasonable  agreement 
with  the  values  by  Gali  et  al.  [14]  (0.50  eV  for  MgAi  and  0.97  eV  for  BeAi)  that  are  also 
obtained  by  HSE  calculations  but  with  a  different  supercell  size.  At  this  point,  both  of 
experimental  and  theoretical  studies  indicated  that  MgAi  and  BeAi  are  not  good  candidates 
for  p- type  doping  in  AIN.  Therefore,  different  approaches  shall  be  taken.  One  of  the  top 
on  the  list  is  by  co-doping  which  Gali  et  al.  already  suggested  a  possible  reduction  in  the 
ionization  energy  with  Mg  and  O  co-doping  that  could  lead  to  MgAi-ON-MgAi  complexes. 


Fig.  2  Formation  energy  as  a  function  of  Fermi  level  for  MgAi  and  BeAi  in  wurtzite 
AIN  under  Al-rich  conditions.  The  Fermi  level  is  shown  in  the  range  from  0  to  5.64; 
corresponding  to  the  HSE  calculated  band  gap  of  AIN.  The  crossing  point  between 
the  0  and  -1  charge  state  defines  the  ionization  energy. 
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Interstitial  aluminum  and  aluminum  vacancy  pairs  (AF-FaQ  in  AIN 

In  the  previous  program,  we  have  studied  the  similar  type  of  pairs  in  GaN,  i.e. 
Garhca-  These  types  of  defect  complexes  are  called  Frenkel  pairs  which  are  believed  to 
fonn  when  the  crystals  are  irradiated  by  electrons  or  particles.  This  is  because,  after 
electron  irradiation  process,  interstitials  and  the  corresponding  vacancies  are  expected. 
Depending  on  the  energy  of  the  irradiated  electron,  A1  or  N  can  be  selectively  knocked 
out  of  their  sites.  Since  Al,  is  a  positive  center  and  Vai  is  a  negative  center,  they  attract 
each  other.  When  they  are  in  an  appropriate  proximity  the  pair  formation  could  be 
expected.  There  are  many  experimental  works  on  the  electronic  and  optical  properties  of 
samples  after  electron  irradiation  in  GaN.  The  attempts  to  detect  these  defects  have  been 
carried  out. [17-21]  In  our  previous  program,  we  studied  the  Ga,  -  VGa  pairs  in  GaN.  We 
found  that  Ga,  can  occur  in  3+,  2+,  and  1+  charge  states,  depending  on  the  Fermi  energy 
of  the  sample.  F(ia  can  occur  in  3-,  2-,  1-,  and  neutral  charge  states.  Because  Ga,-  and 
V Ga  have  opposite  charge,  they  have  Coulom  attraction  to  each  other.  However,  if  the 
two  defects  are  placed  next  to  each  other,  they  are  spontaneously  recombined.  Even  if 
Ga,  is  placed  at  the  next  nearest  neighbor  of  F(ia,  Ga,  can  still  spontaneously  diffuse  via 
the  kick-out  mechanism  and  the  defects  are  annihilated.  The  first  stable 
Ga,  -  VGa  configuration  occurs  when  the  pair  is  separated  by  ~  4  A  along  the  hexagonal 
channel  in  which  Ga,  is  placed  at  octahedral  site.  This  configuration  was  called 
“1st  Ga,  -VGa  ”  in  the  previous  program  and  will  be  called  Ga-FP(l)  here.  The  next  stable 

configuration  with  the  higher  energy  [Ga-FP(2)]  is  found  when  the  Ga,  is  placed  at  the 
next  nearest  octahedral  site  along  the  hexagonal  channel  with  a  distance  of  6.4  A  away 
from  the  ideal  Foa  position.  These  two  pairs  have  the  binding  energy  (compared  the 
energy  of  the  pair  to  the  separated  Ga,-  and  Foa)  of  3.3  and  2.4  eV,  respectively.  We 
found  that  for  further  separation,  the  binding  energy  is  immediately  reduced  to  the  order 
of  0.1  eV. 

In  this  program,  we  performed  similar  study  in  the  case  of  AIN.  Electron 
irradiation  experiments  are  one  of  the  important  means  to  intentionally  create  defects  in 
nitrides  and  oxides  [17,  22].  Since  the  created  defects  can  be  both  interstitials  and 
vacancies,  i.e.,  Al„  N,-,  Vai  and  Fn.  In  this  program,  we  study  the  Al-Frenkel  pairs  (Al- 
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FPs)  which  consist  of  Al,  and  Vai.  We  applied  the  density  functional  theory  (DFT),  as 
implemented  in  the  VASP  code  [1]  with  local  density  approximation  (LDA)  in  our 
calculations.  The  interaction  between  ions  and  electrons  is  described  using  the  projector 
augmented  wave  method  [23].  The  cutoff  energy  for  plane-wave  expansion  of  400  eV  is 
used.  The  96-atom  supercell  is  employed  for  simulation  of  defects  and  the  shifted  2x2x2 
grid  of  Monk-FIorst  Pack  special  k  point  is  used  for  integrations  over  the  Brillouin  zone. 
All  atoms  in  supercell  are  relaxed  until  the  residual  force  is  less  than  0.01  eV/A.  All 
calculations  were  performed  at  the  calculated  lattice  constant  of  bulk  wurtzite  AIN  of  a  = 
3.09  A  and  da  =  1.600. 


Fig.  3  Local  atomic  geometry  of  (a)  the  closest  Al-Frenkel  pair  [Al-FP(l)]  and  (b) 
the  further  separation  Al-Frenkel  pair  [A1-FP(2)]. 

Figure  3  (a)  shows  the  local  atomic  geometry  of  the  first  bound  Al-Frenkel  pair, 
Al-FP(l)  and  (b)  the  second  bound  Al-Frenkel  pair,  A1-FP(2)  which  is  slightly  higher  in 
energy.  The  configurations  of  these  two  Al-Frenkel  pairs  are  very  similar  to  the 
corresponding  Ga-Frenkel  pairs  in  GaN.  Similar  to  the  case  of  Ga-FP,  here,  we  also 
observed  large  relaxations  of  Al  atoms  that  involve  in  Frenkel  pair  stabilization.  The 
formation  energies  for  separated  Al,-  and  Vai  as  well  as  Al-FP(l)  and  A1-FP(2)  are  shown 
in  Fig.  4.  From  the  formation  energy  plot,  we  can  see  that  Al,  is  a  shallow  triple  donor. 
Vai  is  a  deep  triple  acceptor  with  the  transition  levels  s(0/-),  s(-/2-)  and  s(2— /3— )  at  0.46 
eV,  0.82  eV,  and  1.22  eV,  respectively.  The  Al-Frenkel  pairs  can  occur  in  1+,  neutral 
and  2-  charge  states.  For  Al-FP(l),  the  transition  levels  s(+/0)  and  e(0/2— )  are  at  0.53  eV 
and  4.08  eV,  respectively.  For  A1-FP(2),  the  transition  levels  s(+/0)  and  s(0/2— )  are 
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slightly  higher  at  0.56  eV  and  4.28  eV,  respectively.  Table  II  shows  the  binding  energy 
of  Al-Frenkel  pairs  in  neutral  charge  state.  These  binding  energies  are  comparable  to  the 
case  of  Ga-FP  in  GaN.  The  recombination  barrier  of  Al-FP(l)  pair  in  the  neutral  charge 
state  is  ~0.2  eV  which  is  similar  to  the  recombination  barrier  for  Ga-Frenkel  pair  in  GaN. 
This  suggested  that  the  pair  is  unstable  and  would  recombine  at  rather  low  temperatures. 


Fig.  4  Formation  energies  as  a  function  of  Fermi  level  for  Al„  Vai,  and  Al-Frenkel 
pairs;  Al-FP(l)  and  A1-FP(2)  under  Al-rich  conditions.  The  zero  of  Fermi  level 
corresponds  to  the  top  of  the  valence  band.  The  dash  line  is  the  sum  of  formation 
energies  of  the  separated  Al,  and  V\\. 

Table  II  Calculated  binding  energy  of  Al-Frenkel  pairs  (Al-FPs),  initial  distance 
(d i„t)  and  relaxed  distance  (<7rei)  between  Al,  and  VM  in  Al-FP(l)  and  A1-FP(2) 


£bmd  (eV) 

dint  (-A.) 

dre\  (^) 

Al-FP(l) 

2.63 

3.91 

3.78 

A1-FP(2) 

1.71 

6.26 

6.16 
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Resolving  the  coordination  and  neighboring  species  of  In  in  InN  and  In?Cfi  alloys 

We  have  previously  studied  the  local  structures  of  indium  oxynitride  (InOxNi.x) 
with  varied  O  contents  have  been  studied  by  using  the  combination  of  first-principles 
calculations  and  synchrotron  x-ray  absorption  near  edge  structures  (XANES)  [24].  The 
InOxN].x  samples  were  characterized  by  indium  /.3-edgc  XANES  measurements  at 
Synchrotron  Light  Research  Institute  in  Thailand.  The  measurements  were  compared 
favorably  with  the  simulations  of  wurtzite  InN,  wurtzite  InN0.60o.4,  and  bixbyite  ImCf 
crystals.  Because  the  spectrum  of  the  sample  containing  -40%  of  O  (InNo.60o.4)  is  almost 
identical  to  that  of  pure  InN,  we  concluded  that  O  can  substitute  on  the  N  site  in  the 
wurtzite  structure  up  to  40%  before  the  crystal  structure  changed. 

In  this  program,  we  focus  our  attentions  on  understanding  the  source  of  the 
features  in  the  calculated  XANES.  In  the  calculations  we  have  full  control  to  place  the 
neighboring  atom  around  the  In  and  can  simulate  the  spectrum  for  each  case.  For 
wurtzite  InN,  each  In  atom  is  surrounded  by  four  N.  When  O  substitute  for  N  in  the 

wurtzite  structure,  an  In  atom  can  have  one,  two,  three  or  four  O  substitute  for  its  N 

neighbors.  We  have  simulated  for  each  specific  case  and  found  that  the  general  features 
of  the  spectrum  remain  similar  to  that  of  pure  InN  except  that  the  first  peak  becomes 
lower  as  more  O  substitute  for  N.  For  ImCb  which  is  a  Bixbyite  structure,  In  atoms  are 

surrounded  by  six  O  neighbors.  We  found  that  the  XANES  spectrum  of  In  in  the 

Bixbyite  structure  is  clearly  differed  from  that  of  wurtzite.  We  have  also  studied  the 
effects  of  replacing  the  O  neighbor  by  N.  To  understand  the  source  of  the  changes  in  the 
XANES  spectra  in  each  case,  the  decomposed  electron  density  of  states  are  studied. 
Further  details  of  this  work  can  be  found  in  an  attached  manuscript  [P3] . 
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Fig.  5  The  local  structures  around  In 
atoms  in  (a)  bixbyite  ln203,  (b) 
wurtzite  InN,  (c)  wurtzite  InN  with 
four  surrounding  N  atoms  replaced  by 
O  atoms,  and  (d)  wurtzite  InN  with  one 
surrounding  N  atom  replaced  by  an  O 
atom  (left)  and  bixbyite  ln203  with  one 
O  atom  replaced  by  a  N  atom  (right). 
All  bond  distances  are  given  as  a 
percentage  difference  from  an  average 
ln203  bond  distance  (c/ca|c  =  2.170  A). 


Photon  Energy  (eV)  Photon  Energy  (eV) 


Fig.  6  (a)  The  simulated  In  L3-edge  XANES  spectra  of  In  atoms  in  wurtzite  InN 
(bottom),  in  wurtzite  InN  with  one  N  replaced  by  O  (middled),  and  in  wurtzite  InN 
with  four  N  replaced  by  O  (top),  (b)  The  spectra  In  atoms  in  bixbyite  ln203  (bottom 
three  curves)  and  in  bixbyite  ln203  with  one  O  replaced  by  N  (top). 
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Hydrogen  in  ImCf 

In  addition  to  the  study  of  InN  and  In2C>3  alloys,  we  have  also  studied  H  in  In2C>3  with 
the  final  goal  to  study  H  in  the  alloys  (H  in  InN  have  been  previously  studied  [25]). 
In2C>3  is  one  of  a  few  transparent  materials  that  can  be  doped  to  a  very  high  carrier 
concentration  without  significantly  degrading  their  transparencies.  The  most  widely  used 
transparent  conductor  is  tin-doped  In203,  known  as  Indium  Tin  Oxide  (ITO),  which  is 
used  as  a  transparent  electrode  for  optoelectronic  devices  such  as  solar  cells  and  light 
emitting  diodes  (LEDs).  Recently,  Koida  et  al.  showed  that  hydrogen-doped  In203  can 
have  mobility  exceeding  100  cm2/Vs  at  carrier  densities  of  ~1020  cm'3  and  exhibits  good 
transparency  even  at  near  infrared  (NIR)  wavelengths  [26].  In  most  semiconductors  H,  is 
amphoteric,  acting  as  a  donor  in  p- type  samples  and  as  an  acceptor  in  //-type  samples.  If 
H,  is  amphoteric,  it  can  never  be  the  cause  of  conductivity  since  it  self-compensates. 
There  are  exceptions,  though.  For  example,  in  ZnO  H,  was  proposed,  based  on  first- 
principles  calculations  [27],  to  act  exclusively  as  a  donor,  a  prediction  which  was  soon 
confirmed  by  experiments  [28].  Whether  H  is  amphoteric  impurity  like  most 
semiconductors  or  is  acting  exclusively  as  a  donor  in  In2C>3  like  ZnO,  is  detrimental  for 
the  understanding  its  conductivity. 


Fig.  7  Schematic  illustration  of  hydrogen  defects  in  ImCfi,  (a)  ff  at  the  oxygen 
anti-bonding  site,  (b)  H7  at  the  16c  site,  and  (c)  H() . 

In  this  program,  we  employed  first-principles  calculations  using  density  functional 
theory  within  the  local  density  approximation  (LDA),  as  implemented  in  the  Vienna  ab 
initio  simulation  (VASP)  package  [29]  to  study  H  in  In203.  Ion-electron  interactions 
were  treated  using  projector  augmented  wave  potentials  (PAW)  [30].  The  electron  wave 
functions  were  described  using  a  plane-wave  basis  set  with  an  energy  cutoff  of  400  eV. 
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To  study  H  impurities,  a  supercell  approach  with  the  supercell  size  of  80  atoms  was  used. 
We  found  that  H,  is  a  shallow  donor  in  In203,  and  can  therefore  also  be  the  cause  of  77- 
type  conductivity  (similar  to  ZnO).  In  addition,  H  atoms  can  also  occupy  substitutional 
oxygen  sites  (Ho)  and  also  act  exclusively  as  donors  in  this  configuration.  Substitutional 
hydrogen  was  recently  described  as  a  multicenter  bond  and  found  to  be  stable  in  a 
number  of  oxides  and  some  nitride  materials  [31-33].  In  order  to  explore  the  stability  of 
Ho,  we  also  performed  calculations  for  oxygen  vacancies  (Fo)  in  ImCfi.  Our  findings 
confirm  that  hydrogen  acts  as  source  of  /?-type  conductivity  in  ImO-,. 

Further  details  of  this  work  can  be  found  in  an  attached  manuscript  [P6] . 


Fig.  8  Calculated  formation  energies  of  hydrogen  defects  and  oxygen 
vacancies  in  Im03  as  a  function  of  Fermi  energy.  In-rich  conditions  are 
assumed. 
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High  pressure  study  of  GaN,  InN,  and  other  wurtzite  semiconductors 

In  actual  devices,  several  layers  of  different  semiconductors  are  generally 
required  to  be  stacked  up  together.  This  often  leads  to  large  strain  while  most  of  the  first 
principles  studies  are  assumed  strain  free.  In  order  to  pave  the  way  to  include  strain 
effects,  we  study  the  elastic  constants  and  sound  velocities  as  function  of  direction  and 
pressure  in  various  wurzite  crystals.  These  include  GaN  and  InN  as  well  as  other 
wurtzite  semiconductors,  i.e.  SiC,  ZnO,  and  CdSe.  We  did  not  limit  our  study  to  just  III- 
V  compounds  because  we  would  like  to  systematically  study  the  trend  of  the  elastic 
constants  and  sound  velocities  with  respect  to  the  ionicity  of  material. 

In  this  program,  to  obtain  all  total  energies  required  for  the  calculation  of  the 
elastic  constants  we  use  the  Kohn-Sham  density  functional  theory  [34,  35]  in  the  local 
density  approximation  with  the  Perdew-Zunger  parametrization  of  exchange  and 
correlation.  [3  6]  The  total  energies  are  calculated  using  the  fullpotential  linearized  muffin- 
tin  orbital  method  as  described  by  Methfessel  and  van  Schilfgaarde.[37]  Well  converged 
double  basis  sets  are  utilized  and  the  Brillouin  zone  integrations  are  done  using  a  4x4x4 
Monkhorst-Pack  sampling  set.  [3 8] 

The  elastic  constants  and  sound  velocities  in  wurtzite  GaN  and  InN  as  well  as 
other  wurtzite  semiconductors  (SiC,  ZnO,  and  CdSe)  as  functions  of  pressure  are 
calculated.  We  found  interesting  nonlinear  behavior  of  the  elastic  constants  and  sound 
velocities  with  pressure  that  is  quite  different  in  the  different  materials.  We  also  showed 
that  the  elastic  constants  are  lower  in  materials  with  higher  ionicity.  Generally,  the  sound 
velocities  follow  the  same  trend  as  the  elastic  constants,  i.e.,  lower  in  materials  with 
higher  ionicity.  The  pressures  at  which  the  C44  and  C66  shear  elastic  constants  cross  are 
investigated.  When  becomes  lower,  it  should  be  easier  to  excite  the  symmetry 
breaking  strain  mode  related  to  the  wurtzite-to-rocksalt  phase  transition.  We  found  that 
the  crossover  pressure  is  generally  higher  than  the  experimental  transition  pressure  and 
explained  it  by  the  fact  that  finite  temperature  fluctuations  make  it  possible  to  excite  the 
required  modes  even  before  the  two  sound  velocities  cross.  We  proposed  the  average  of 
the  crossover  and  calculated  equilibrium  transition  pressures  to  be  a  good  estimate  for  the 
experimental  transition  pressure  in  cases  where  there  is  a  significant  transition  enthalpy 
barrier.  Of  course,  the  details  near  the  transition  point,  the  height  of  the  barrier,  and  the 
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coupling  of  the  in  and  out  of  plane  acoustic  modes,  which  are  both  required  to  reach  the 
transition  point,  are  important  factors  for  the  crystal  phase  transformation. 

Further  details  of  this  work  can  be  found  in  an  attached  manuscript  [P2]. 


Pressure  (GPa) 


Fig.  9  Calculated  pressure  dependence  of  the  elastic  constants  in  (a)  SiC,  (b) 
GaN,  (c)  InN,  (d)  ZnO,  and  (e)  CdSe.  In  (b)  and  (c),  the  square  symbols  show 
the  calculated  results  by  Lepkowski  et  al.  [39] 
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Properties  of  phase  separated  Zn3P?  in  ZnO 

Acceptor  doping  of  ZnO  has  proved  highly  challenging.  Various  dopants  have 
been  tried  with  limited  experimental  reports  of  success.  While  p- type  conductivity  has 
been  reported  in  ZnO  doped  with  group-V  acceptors  such  as  P,  As  and  Sb,  computational 
studies  have  systematically  produced  large  values  for  the  corresponding  acceptor 
ionization  energies,  rendering  it  unlikely  that  these  elements  would  act  as  shallow 
dopants.  In  this  program,  propose  a  mechanism  that  may  give  rise  to  the  measurements 
of  hole  conductivity  in  acceptor-doped  samples.  When  P  is  used  as  an  acceptor  dopant, 
the  intent  is  to  introduce  large  concentrations  of  the  impurity  into  ZnO.  Under  such 
conditions,  the  solubility  limit  may  be  exceeded  and  a  phase  separated  Zn3P2  might  be 
fonned. 

Based  on  on  density  functional  theory  in  the  local  density  approximation,  we 
investigated  phase  separated  Zn3P2  in  ZnO.  (Ultrasoft  pseudopotentials  with  a  plane- 
wave  basis  set  -  energy  cutoff  300  eV  -  as  implemented  in  the  VASP  code  [29]  were 
used.)  We  found  that  the  presence  of  Zn3P2  can  lead  to  measured  hole  conduction  in  two 
ways.  First,  the  Zn3P2  itself  exhibits  a  tendency  towards  p- type  conductivity.  Second, 
we  demonstrated  that  at  the  interface  between  ZnO  and  the  Zn3P2  there  is  a  strong 
tendency  for  a  hole  accumulation  layer  to  fonn. 

Further  details  of  this  work  can  be  found  in  an  attached  manuscript  [PI]. 
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Oxygen  defects  in  CdTe 

Recently,  local  vibrational  inodes  (LVM)  related  to  oxygen  impurities  in  CdTe 
have  been  experimentally  studied  using  Fourier  transform  infrared  (FTIR)  spectroscopy 
[40,  41].  Depending  on  the  growth  and  doping  conditions,  low  temperature  FTIR 
measurements  reveal  three  LVMs  related  to  oxygen:  one  low  frequency  mode  at  350  cm'1 
and  two  high  frequency  modes  at  V\  =  1097  and  vy  =  1108  cm'1.  As  the  measuring 
temperature  is  increased,  the  two  high  frequency  modes  merge  into  one  at  room 
temperature.  Chen  et  al.  assigned  the  low  frequency  mode  to  an  oxygen  substituting  on 
the  Te  site  (0Te),  and  the  two  high  frequency  modes  to  a  complex  formed  by  an  oxygen 
substituting  on  a  Te  site  and  a  neighboring  Cd  vacancy  (Oie-Fcd)  [40,  41]. 
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Fig.  11  The  local  structures  of  (a)  Bulk  CdTe,  (b)  0Te,  and  (c)  0Te-Fcd 
complex  in  CdTe. 

To  understand  these  vibrational  modes,  we  have  investigated  O-related  centers  in 
CdTe  using  density  functional  calculations.  Our  work  is  based  on  density  functional 
theory  within  the  local  density  approximations  (LDA).  For  the  electron-ion  iteractions, 
we  used  projector  augmented  wave  potentials  [23],  as  implemented  in  the  VASP  code 
[42,  43],  with  an  energy  cutoff  of  500  eV  for  the  plane- wave  basis  set.  The  LVM 
frequencies  of  Oie  and  Oie-Fcd  defects  were  calculated  by  using  force-response  matrix 
(dynamical  matrix)  calculations.  The  harmonic  approximation  is  used.  We  found  that  the 
calculated  LVM  of  0|e  is  consistent  with  the  experimental  value.  Flowever,  the 
calculated  LVM  frequencies  of  0Te-  Fed  are  much  lower  than  the  V\  and  Vj  modes 
observed  by  Chen  et  al.  The  calculations  clearly  showed  that  the  0Te-Fcd  complex  is  not 
the  correct  model  to  explain  the  observed  LVMs.  Preliminary  results  of  our  work  have 
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been  published  as  a  Comment. [44]  In  this  program,  we  laid  out  the  computational 
method  and  details  of  the  LVMs.  The  two  high  frequencies  observed  by  Chen  et  al. 
should  be  attributed  to  other  defects  related  to  the  introduction  of  O  in  CdTe,  not  to  0Te- 
VCd.  We  note  that  the  high  value  of  the  frequency  makes  hydrogen  (which  could  be 
related  to  the  introduction  of  O  during  growth)  a  prime  candidate. 

Further  details  of  this  work  can  be  found  in  an  attached  manuscript  [P4]. 


Frequency  (cm  ) 


Fig.  12  The  sum-square  of  the  three  eigenvector  components  associated  with 
the  O  atom  for  all  vibrational  modes  for  (a)  Oie,  and  (b)  the  Oie-Fcd  complex 
(in  a  64-atom  CdTe  supercell).  While  the  calculated  frequency  of  Oie  is  in 
agreement  with  the  observed  value,  the  calculated  frequencies  of  Oie-Fcd  are 
much  lower  than  the  observed  values  of  -1100  cm'1;  suggesting  that  the 
model  proposed  by  Chen  et  al.  is  invalid. 
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Effects  of  unit  cell  size  for  studying  alloy  properties 

Studying  alloys  by  first  principles  calculations  is  a  very  challenging  task.  This  is 
because,  for  a  given  alloy  concentration,  there  are  virtually  infinite  ways  to  arrange 
atoms.  Since  first  principles  calculations  can  handle  a  system  of  the  order  of  100  atoms 
and  not  much  more,  alloys  are  generally  studied  by  the  supercell  approach.  This  means 
the  supercell  unit  is  repeated  and,  strictly  speaking,  the  system  is  ordered.  In  order  to 
simulate  disordered  alloys,  many  configurations  have  to  be  calculated  and  some 
averaging  schemes  have  to  be  employed. 
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Fig.  13  The  contour  plots  of  the  calculated  bandgap  ranges  of  pseudoternary 
AlGalnP  as  a  function  of  alloy  compositions  (u  =  x  +  y/ 2  and  v  =  V3y/2). 

Each  figure  shows  the  results  of  the  calculation  with  a  different  Amax. 

One  of  the  most  studied  properties  of  alloys  is  the  band  gap.  For  III-V  alloys,  the 
cation  alloys,  for  examples  AlxGayIni_x-yN,  AlxGayIni_x_yAs,  or  AlxGayIni_x_yP  are  widely 
studied  both  experimentally  and  theoretically.  It  is  interesting  that  experimental  results 
from  different  growth  techniques  are  often  in  disagreement  even  when  the  same  or 
similar  alloy  compositions  are  compared.  We  have  suggested  that  one  of  the  cause  of 
these  disagreements  might  be  the  fact  that  samples  from  different  techniques  have 
different  alloy  orderings  [45].  In  this  program,  we  focused  our  attentions  in  studying  a 
complete  set  of  theoretical  possible  arrangements  within  a  given  cell  size.  These 
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complete  results  allow  us  to  investigate  how  the  band  properties  relate  to  the  atomic 
arrangements.  Note,  however,  that  many  of  the  configurations  studied  might  have  very 
high  energy  and  difficult  to  fabricate.  While  a  larger  cell  size  provides  more  possible 
configurations  -  allowing  a  more  complete  study  -  and  finer  composition  steps,  the 
computational  demand  is  increased  rapidly.  Therefore,  it  is  important  to  investigate  the 
size  of  unit  cell  that  sufficiently  gives  the  converged  results.  In  this  work,  we  used  the 
quaternary  alloys  AlxGayIni.x.yP  as  an  illustration  case.  The  general  results  should, 
however,  be  applied  to  other  alloy  systems.  We  found  that  it  is  sufficient  to  use  the 
supercell  size  as  small  as  5  cation  atoms  to  yield  reasonable  results;  providing  that  all 
possible  cell  shapes  are  included.  The  bandgap  range  of  AlGalnP  alloy  is  reported  (in  the 
fonn  of  a  polynomial  function).  It  is  found  that  the  bandgap  range  for  AlGalnP  alloy 
system  is  large  (~0.8  eV)  for  the  compositions  near  the  middle  of  the  composition  space 
(slightly  shifted  toward  the  low  A1  composition  side). 

Further  details  of  this  work  can  be  found  in  an  attached  manuscript  [P5] . 
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We  report  first-principles  calculations  and  interface  simulations  for  Zn3P2,  a  compound  that  may 
form  during  doping  of  ZnO  with  phosphorous.  While  P  is  a  deep  acceptor  in  ZnO  and  thus  unable 
to  produce  p-type  conductivity,  we  show  that  hole  accumulation  can  occur  at  ZnO/Zn3P2  interfaces 
due  to  the  unusual  valence-band  alignment  between  the  two  materials.  This  provides  an  explanation 
for  the  hole  conductivity  that  has  been  observed  in  Hall  measurements  on  phosphorous-doped 
ZnO.  ©  2010  American  Institute  of  Physics.  [doi:10. 1063/1. 3481069] 


ZnO  is  a  highly  promising  material  for  a  wide  variety  of 
applications,  including  window  layers  for  solar  cells,  trans¬ 
parent  transistors,  and  sensors.  Based  on  its  band  gap  (3.4 
eV)  and  large  exciton  binding  energy  (60  meV)  it  would  also 
be  an  excellent  material  for  light  emitters.  Acceptor  doping 
of  ZnO  has  proved  highly  challenging.  There  are  by  now  a 
large  number  of  reports  in  the  literature  (too  numerous  to 
cite)  claiming  observations  of  p-type  conductivity  in  ZnO 
doped  with  different  elements.  The  lack  of  follow-up  reports 
and  the  scarcity  of  device  demonstrations  based  on  p-n  junc¬ 
tions  indicate  that  there  are  serious  problems  with  reliability 
and  reproducibility.  In  addition,  fundamental  questions  have 
arisen  about  the  ability  of  some  of  the  acceptor  dopants  to 
actually  yield  hole  conduction.  In  particular,  while  p-type 
conductivity  has  been  reported  in  ZnO  doped  with  group-V 
acceptors  such  as  P,1-4  As,5  and  Sb,6  computational  studies7 
have  systematically  produced  large  values  for  the  corre¬ 
sponding  acceptor  ionization  energies,  rendering  it  unlikely 
that  these  elements  would  act  as  shallow  dopants.  In  the  case 
of  nitrogen,  estimates  of  ionization  energies  have  been  in¬ 
conclusive,  although  recent  computations  indicate  a  very 

g 

deep  level.  Nitrogen  is  the  element  for  which  the  most  fa¬ 
vorable  results  have  been  obtained;7  11  still,  the  questions 
about  reproducibility  and  reliability  also  apply  in  this  case. 

If  these  elements,  when  incorporated  into  ZnO,  do  not 
actually  act  as  shallow  acceptors,  then  why  do  Hall  measure¬ 
ments  yield  a  signal  that  seems  indicative  of  hole  conductiv¬ 
ity?  One  possible  cause,  recently  discussed  by  Ohgaki 
et  al.  and  by  Bierwagen  et  al.  is  the  presence  of  inhomo¬ 
geneities  in  the  sample.  It  was  demonstrated,  based  both  on 
modeling  and  on  explicit  measurements  for  a  known  n-type 
sample,  that  lateral  inhomogeneities  in  carrier  concentrations 
can  result  in  an  incorrect  assignment  of  the  carrier  type. 
Other  mechanisms  for  explaining  the  hole  conductivity  have 
been  proposed1417  but  have  no  direct  experimental  confirma¬ 
tion. 

In  this  letter,  we  propose  and  document  another  mecha¬ 
nism  that  may  give  rise  to  spurious  measurements  of  hole 
conductivity  in  acceptor-doped  samples.  The  mechanism  will 
be  illustrated  with  the  example  of  P  but  it  may  also  occur  in 
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the  case  of  the  other  group-V  impurities  (N,  As,  and  Sb). 
When  P  is  used  as  an  acceptor  dopant,  the  intent  is  to  intro¬ 
duce  large  concentrations  of  the  impurity  into  ZnO.  Under 
those  conditions,  the  solubility  limit  may  be  exceeded,  an 
issue  that  has  been  extensively  discussed  in  the  context  of 
doping  of  other  semiconductors  (see  Refs.  16  and  17  for  a 
discussion  of  acceptors  in  ZnSe  and  ZnTe  and  Refs.  18  and 
19  for  GaN).  Solubility  limits  arise  because  it  becomes  more 
favorable  for  the  impurity  to  incorporate  into  an  alternate 
phase  rather  than  on  a  substitutional  lattice  site  inside  the 
semiconductor.  The  competing  phase  would  typically  mani¬ 
fest  itself  as  precipitates  inside  the  host  material.  Candidates 
for  alternate  phases  can  be  identified  by  examining  which 
compounds  can  be  formed  between  the  impurity  and  the  host 
element.  In  the  case  of  group-V  impurities  in  ZnO,  the  prime 
candidates  are  the  Zn3Y2  compounds  (where  Y=N,  P,  As,  or 
Sb).  We  therefore  perform  an  investigation  of  such  com¬ 
pounds  and  their  impact  when  they  are  incorporated  as  pre¬ 
cipitates  in  ZnO.  We  will  see  that  the  presence  of  such  com¬ 
pounds  can  lead  to  measured  hole  conduction  in  two  ways. 
First,  the  compounds  themselves  exhibit  a  tendency  toward 
p-type  conductivity.  Second,  we  will  demonstrate  that  at  the 
interface  between  ZnO  and  the  Zn3Y2  compound  there  is  a 
strong  tendency  for  a  hole  accumulation  layer  to  form. 

Our  investigations  of  ZnO  and  Zn3P2  are  based  on  first- 
principles  calculations  using  density  functional  theory  in  the 
local  density  approximation.  Ultrasoft  pseudopotentials  with 
a  plane-wave  basis  set  (energy  cutoff  300  eV)  as  imple¬ 
mented  in  the  VASP  code20  were  used.  Bulk  calculations  were 
performed  for  Zn3P2  in  the  PAH  nmc  structure  with  a  40- 
atom  unit  cell.  We  obtained  a  calculated  lattice  constant  of 
7.897  A,  which  compares  favorably  with  the  experimental 
value21  of  8.097  A.  A  2X2X2  Monkhorst-Pack  k-point 
sampling  was  used  in  the  Brillouin-zone  integration. 

Insight  into  band  alignments  as  well  as  doping  tenden¬ 
cies  of  the  material  can  be  provided  by  investigating  the 
behavior  of  interstitial  hydrogen.-2  We  performed  supercell 
calculations  for  a  single  hydrogen  interstitial  in  Zn3P2  using 
40-atom  supercells  and  2X2X2  k-point  sampling.  The  for¬ 
mation  energy  of  H?  ( q=- ,  0,  or  +)  is  calculated  according  to 
the  definition  given  in  Ref.  23. 

Consistent  with  general  tendencies,--  H  in  the  positive 
charge  state  prefers  a  lattice  location  in  a  P-antibonding  site. 
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FIG.  1 .  (Color  online)  Calculated  formation  energies  of  hydrogen  intersti¬ 
tials  in  Zn3P2  as  a  function  of  Fermi  level  (referenced  to  the  VBM).  H2 
molecules  are  assumed  as  the  reservoir. 


while  in  the  negative  charge  state  it  prefers  a  Zn-antibonding 
site.  The  neutral  charge  state  is  never  stable,  as  can  be  seen 
in  Fig.  1.  The  +/-  transition  level  is  found  at  0.12  eV  above 
the  valence-band  maximum  (VBM). 

Assuming  the  universal  alignment  discussed  in  Ref.  22 
is  valid,  the  position  of  the  VBM  with  respect  to  the  vacuum 
level  is  expected  to  be  at  approximately  -4.50-0.12 
=  -4.62  eV.  This  is  quite  high,  on  an  absolute  energy  scale, 
and  will  have  consequences  for  the  electronic  behavior  of 
interfaces  between  Zn3P2  and  ZnO,  as  discussed  below.  Ma¬ 
terials  in  which  the  VBM  occurs  at  a  high  energy  with  re¬ 
spect  to  vacuum  generally  exhibit  a  tendency  toward  p-type 
conductivity,  i.e.,  impurities  that  act  as  shallow  acceptors 
will  be  more  easily  incorporated,  and  compensation  effects 
due  to  donors  suppressed.  This  tendency  toward  /j-type  con¬ 
ductivity  may,  in  itself,  explain  why  P-doped  ZnO  in  which 
Zn3P2  precipitates  occur  could  exhibit  measurable  hole  con¬ 
duction.  However,  as  we  will  discuss  below,  it  is  not  actually 
necessary  for  the  Zn3P2  itself  to  be  p  type. 

Another  means  of  producing  insight  into  the  band  align¬ 
ment  is  by  performing  calculations  for  a  Zn3P2  surface, 
which  produce  the  position  of  the  VBM  with  respect  to  the 
vacuum  level.  In  principle,  such  alignments  are  sensitive  to 
the  details  of  the  surface  but  for  the  qualitative  purposes  of 
the  present  work  this  is  not  a  limitation.  We  performed  sur¬ 
face  calculations  in  a  repeated-slab  geometry  with  a  slab 
thickness  of  nine  layers  (approximately  16  A)  along  the  non¬ 
polar  [010]  surface  and  a  vacuum  spacing  of  12  A  in  be¬ 
tween.  On  the  unrelaxed  surface,  the  VBM  is  found  at  —4.38 
eV  with  respect  to  the  vacuum  level.  This  value  changes  to 
—4.48  eV  when  relaxations  are  included,  indicating  the  value 
is  not  very  sensitive  to  such  effects.  Gratifyingly,  the  latter 
value  is  in  good  agreement  with  the  value  for  the  Zn3P3 
VBM  position,  at  —4.62  eV,  obtained  from  the  hydrogen- 
level  alignment. 

We  now  turn  to  simulations  of  the  ZnO/Zn3P2  interface. 
The  purpose  of  this  simulation  is  not  to  provide  an  atomistic 
depiction  of  the  interface;  indeed,  because  of  the  sizeable 
lattice  mismatch,  creating  a  coherent  interface  between  ZnO 
and  Zn3P2  would  be  challenging.  As  we  will  see,  however, 
our  conclusions  do  not  depend  on  details  of  interfacial  struc¬ 
ture  but  only  on  qualitative  values  of  band  alignment.  There¬ 
fore  use  of  “natural  band  alignments,”  in  this  case  based  on 
the  estimates  obtained  above  for  hydrogen-level  alignments 


FIG.  2.  (Color  online)  (Bottom)  Band  diagram  for  a  ZnO/Zn3P2  interface. 
(Top)  Corresponding  density  of  holes. 


and  surface  calculations,  will  suffice.  We  focus  on  a  simula¬ 
tion  of  band  positions  in  the  vicinity  of  the  interface. 

Our  estimates  above  placed  the  VBM  of  Zn3P2  at  around 
4. 5-4. 6  eV  below  the  vacuum  level.  Combined  with  infor- 
mation  about  the  VBM  of  ZnO,  we  estimate  a  valence- 
band  offset  between  the  two  materials  of  about  3.4  eV. 
Interface-specific  effects  may  affect  this  value  at  an  actual 
heterojunction;  however,  we  will  demonstrate  that  our  con¬ 
clusions  are  not  sensitive  to  the  precise  value  of  this  offset 
(with  a  very  wide  tolerance).  We  note  that,  with  a  ZnO  band 
gap  of  3.40  eV,  this  estimated  VB  offset  causes  the  VBM  of 
Zn3P2  to  lie  at  approximately  the  same  level  as  the 
conduction-band  minimum  (CBM)  of  ZnO.  This,  in  its  own 
right,  indicates  that  at  an  interface  there  may  be  a  strong 
tendency  for  electrons  from  the  Zn3P2  VBM  to  transfer  to  the 
ZnO  CBM,  thus  leading  to  the  formation  of  holes  in  Zn3P2. 

We  have  backed  up  these  qualitative  insights  with  ex¬ 
plicit  Schrodinger-Poisson  simulations”4  of  a  ZnO/Zn3P2  in¬ 
terface.  For  the  purpose  of  the  simulations  some  assumptions 
need  to  be  made  about  doping  of  the  two  materials.  (We  will 
show  below  that  our  results  are  not  sensitive  to  the  details  of 
these  assumptions.)  On  the  ZnO  side,  we  assume  that  P  is 
incorporated  as  a  deep  acceptor  in  ZnO  with  an  ionization 
energy  of  0.93  eV  (Ref.  7).  The  simulation  in  Fig.  2  was 
performed  assuming  a  deep  acceptor  concentration  of 
1019  cm-3.  With  regards  to  doping  of  Zn3P2,  as  noted  above 
we  expect  this  material  to  have  a  tendency  to  be  p  type; 
however,  for  the  purposes  of  the  simulation  we  assumed  a 
“worst-case”  scenario,  in  which  Zn3P2  is  n  type  (with  a  do¬ 
nor  concentration  of  1017  cm-3),  thus  making  it  very  clear 
that  the  appearance  of  holes  at  the  interface  does  not  rely  on 
the  p-type  nature  of  Zn3P2.  A  band  gap  of  1.5  eV  was  as¬ 
sumed  for  Zn3P2  (Ref.  25). 

Figure  2(a)  clearly  shows  that  a  hole  accumulation  layer 
is  formed  at  the  interface.  Because  the  Zn3P2  VBM  lies  at  a 
high  energy,  electrons  in  the  Zn3P2  VB  have  a  tendency  to 
transfer  to  the  ZnO,  where  they  can  find  lower-energy  states. 
If  the  ZnO  were  undoped  or  n  type,  these  electrons  would  go 
into  the  ZnO  conduction  band  (CB);  the  interface  would  then 
exhibit  mixed  conduction.  But  if  the  ZnO  is  doped  with  deep 
acceptors,  the  electrons  can  compensate  those  acceptors  and 
free  electrons  are  removed  from  the  system.  Figure  2(b) 
shows  that  a  hole  accumulation  layer  forms  within  a  distance 
of  a  few  nanometer  from  the  interface,  and  with  a  density 
exceeding  102°  cm-3.  The  corresponding  charge  density  of 
the  sheet  is  on  the  order  of  1013  cm-2  (Fig.  3).  We  suggest 
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FIG.  3.  (Color  online)  Dependence  of  integrated  hole  density  in  the  hole 
accumulation  layer  at  a  ZnO/Zn3P2  interface  on  the  value  of  the  band  offset 
between  ZnO  and  Zn3P2.  Zero  offset  (dashed  line)  refers  to  the  situation 
where  the  ZnO  CBM  is  lined  up  with  the  Zn3P2  VBM;  positive  values 
indicate  that  the  VBM  in  Zn3P2  is  raised  with  respect  to  ZnO. 

that  the  presence  of  this  sheet  of  charge,  at  interfaces  embed¬ 
ded  in  the  ZnO  layer,  may  be  responsible  for  measurements 
of  hole  conductivity  in  P-doped  ZnO. 

We  have  performed  extensive  checks  to  ensure  that  our 
qualitative  conclusions  are  not  sensitive  to  the  details  of  our 
simulations.  Decreasing  the  concentration  of  deep  acceptors 
in  ZnO  leads  to  a  decrease  in  the  density  of  the  hole  accu¬ 
mulation  layer,  but  as  noted  above,  a  p-type  layer  continues 
to  form  even  if  the  ZnO  is  undoped  or  n-type.  We  have  also 
checked  the  sensitivity  to  the  band  alignment.  Figure  3,  gen¬ 
erated  using  parameters  identical  to  those  in  Fig.  2  but  with  a 
variable  band  alignment,  clearly  shows  that  a  hole  accumu¬ 
lation  layer  continues  to  form,  even  if  the  Zn3P2  VBM  is 
well  below  the  ZnO  CBM. 

While  our  calculations  are  based  on  P-doped  ZnO,  we 
believe  that  the  conclusions  also  apply  to  ZnO  doped  with 
other  deep  acceptors  such  as  As  or  Sb.  Indeed,  explicit  cal¬ 
culations  for  Zn3As2  showed  that  the  VBM  at  its  (relaxed) 
surface  is  located  at  —4.53  eV,  a  result  very  similar  to  that 
obtained  for  Zn3P2.  Following  the  general  trends  between 
phosphides,  arsenides,  and  antimonides,  we  expect  the 
VBM  of  Zn3Sb2  to  lie  at  least  as  high  as  that  of  Zn3P2  and 
Zn3As2.  The  VBM  of  nitride  compounds  is  expected  to  lie 
lower  than  that  of  P  compounds.  However,  the  fact  that  a 
hole  accumulation  layer  can  form  even  if  the  VBM  of  Zn3P2 
is  2  eV  below  the  CBM  of  ZnO,  as  shown  in  Fig.  3,  is  a 
strong  indicator  that  the  mechanism  described  here  may  also 
apply  in  the  case  of  Zn3N2.  This  could  be  a  potential  expla¬ 
nation  for  the  localized  p-type  conductivity  phenomena  ob¬ 
served  in  Ref.  26. 


In  summary,  we  have  presented  first-principles  calcula¬ 
tions  for  Zn3P2  and  simulations  of  ZnO/Zn3P2  interfaces  that 
provide  a  plausible  explanation  for  the  observed  p- type  con¬ 
ductivity  in  ZnO  doped  with  group-V  acceptors. 
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Elastic  constants  and  sound  velocities  calculated  from  first  principles  as  function  of  pressure  are  presented 
for  wurtzite  SiC,  GaN,  InN,  ZnO,  and  CdSe.  The  Clt  and  C33  elastic  constants,  which  are  involved  in 
longitudinal  sound  waves  along  symmetry  directions,  are  found  to  monotonically  increase  with  pressure.  The 
shear  moduli  C44  and  C66,  which  are  involved  in  transverse  sound  waves  along  symmetry  directions,  either 
decrease  with  increasing  pressure  or  initially  increase  from  zero  pressure  but  then  turn  over  and  start  decreas¬ 
ing.  Of  special  interest  is  the  pressure  at  which  the  C44  and  C 66  elastic  constants  cross.  At  this  pressure,  the 
transverse  acoustic  waves  in  the  basal  plane,  which  are  shown  to  be  closely  related  to  the  symmetry  breaking 
strain  component  that  leads  to  the  phase  transition,  become  easier  to  excite  than  the  ones  with  displacement 
along  the  c  axis.  It  is  found  that  this  crossover  pressure  is  an  upper  limit  to  the  actual  phase  transition  pressure. 

The  average  of  the  calculated  equilibrium  transition  pressure  and  the  crossover  pressure  is  proposed  as  a  good 
estimate  for  the  actual  transition  pressure  in  cases  where  the  transition  is  strongly  kinetically  hindered  by  an 
enthalpy  barrier  between  the  two  phases.  This  occurs  for  SiC  and  GaN  and  is  confirmed  with  literature  data  for 
AIN.  For  the  remaining  materials,  all  these  pressures  are  close  to  each  other.  The  trends  of  the  elastic  constants 
and  sound  velocities  with  the  materials’  Phillips  scale  ionicity  are  also  reported. 

DOI:  10.1103/PhysRevB. 82.035201  PACS  number(s);  61.50.Ks,  62.65. +k 


I.  INTRODUCTION 

Most  group-IV,  III-V,  and  II-VI-semiconductors  are  well 
known  to  undergo  a  high-pressure  phase  transition  from  a 
tetrahedrally  bonded  phase  to  an  octahedrally  bonded  phase. 
Most  prior  work  has  focused  on  determining  the  equilibrium 
transition  pressure,  at  which  the  two  phases  have  the  same 
enthalpy.  Recently,  however,  several  studies  have  started  to 
study  the  mechanism  of  transformation  from  one  phase  to  the 
other.  In  particular,  homogeneous  transition  paths  have  been 
proposed  between  wurtzite  and  rocksalt1  3  and  zincblende 
and  rocksalt. 4-7  In  such  mechanisms,  the  transformation  hap¬ 
pens  by  a  continuous  deformation  from  one  structure  to  the 
other  by  imposed  strains.  At  each  strain,  the  internal  coordi¬ 
nates  (i.e.,  atomic  positions)  are  required  to  adjust  to  mini¬ 
mize  the  energy.  For  example,  for  the  wurtzite  to  rocksalt 
transformation,  the  two  strains  correspond  to  a  da  reduction 
perpendicular  to  the  basal  plane  and  a  strain  in  the  basal 
plane  along  the  [1010]  direction  or  any  direction  which  is 
±60°  or  ±120°  from  it,  which  we  can  call  the  bla  direction. 
On  the  experimental  side,  advanced  diffraction  techniques 
utilizing  high  brilliance  synchrotron  x-rays  open  up  the  op¬ 
portunity  to  study  crystal  structures  under  high  pressure  in 
more  detail. 8-11  Although  the  actual  transition  mechanism  is 
likely  to  be  more  complex  than  a  homogeneous  transforma¬ 
tion,  as  shown  for  example  by  molecular  dynamics  studies,12 
the  homogeneous  transformation  path  provides  a  useful  sim¬ 
plified  picture  which  may  in  reality  happen  locally  in  small 


regions  and  determine  the  nucleation  event  of  the  transition. 
Sarasamak  et  al?  studied  the  enthalpy  landscapes  in  this 
space  of  strains  driving  the  homogeneous  transformation. 
Such  enthalpy  landscapes  can  be  constructed  at  different  hy¬ 
drostatic  pressures.  By  constructing  such  enthalpy  land¬ 
scapes,  one  can  gain  insight  in  the  behavior  of  the  system 
along  various  paths  from  one  phase  to  the  other  within  the 
strain  coordinates.  The  local  curvatures  in  different  direc¬ 
tions  for  example  will  tell  us  in  which  direction  the  system 
will  initially  most  easily  deform.  The  elastic  constants  are 
closely  related  to  the  local  shape  of  this  enthalpy  landscape. 
In  fact,  they  essentially  determine  the  curvatures  near  the 
minima  on  this  energy  landscape.  It  is  thus  of  interest  to 
determine  the  elastic  constants  as  function  of  pressure  to  bet¬ 
ter  understand  the  initial  stages  of  phase  transition. 

The  elastic  constants  are  also  closely  related  to  the  sound 
velocities.  Thus  studying  the  sound  velocities  as  function  of 
pressure  provides  possibly  an  experimental  way  to  determine 
the  pressure  dependence  of  the  elastic  constants.  In  fact,  ho¬ 
mogeneous  strains  are  simply  the  long-wavelength  limit  of 
the  acoustic  phonon  modes.  The  relation  between  sound 
waves  in  specific  directions  and  the  strains  involved  in  the 
zincblende  and  wurtzite  to  rocksalt  transition  were  previ¬ 
ously  pointed  out  by  Prikhodko  et  al. 13 

Here,  we  study  the  elastic  constants  and  sound  velocities 
as  function  of  direction  and  pressure  in  various  wurzite  crys¬ 
tals:  SiC,  GaN,  InN,  ZnO,  and  CdSe.  Our  choice  of  materials 
includes  a  IV-IV,  two  III-V,  and  two  II- VI  compounds.  These 
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materials  also  have  varied  transition  pressures.  As  will  be 
seen,  rather  different  behavior  is  obtained  for  different  mate¬ 
rials.  Because  we  have  systematically  studied  several  mate¬ 
rials,  we  can  investigate  the  trend  of  the  elastic  constants  and 
sound  velocities  with  respect  to  the  ionicity  of  material. 

The  paper  is  organized  as  follows.  In  Sec.  II,  we  give  the 
expressions  for  the  sound  velocities  in  different  direction  as 
function  of  the  elastic  constants.  The  calculation  of  elastic 
constants  as  function  of  pressure  is  also  discussed.  Essen¬ 
tially,  they  are  obtained  by  calculating  the  total  energy  as 
function  of  a  set  of  strains.  Details  of  the  computational 
method  used  to  calculate  the  total  energies  from  first  prin¬ 
ciples  are  given  in  Sec.  III.  The  results  for  the  elastic  con¬ 
stants  and  sound  velocities  in  specific  high-symmetry  direc¬ 
tions  as  function  of  pressure  are  presented  in  Sec.  IV.  In  this 
section,  the  trend  of  the  elastic  constants  and  sound  veloci¬ 
ties  with  respect  to  the  ionicity  of  material  with  also  be  pre¬ 
sented.  A  summary  of  the  main  results  is  given  in  Sec.  V. 
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H.  THEORY 
A.  Sound  velocities 

The  acoustic  vibrational  modes  of  a  solid  in  the  long- 
wavelength  limit  are  obtained  from  the  Christoffel 
equation,14 

p(o2Uj  =  Muuh  (1) 

where  p  is  the  mass  density,  u>  the  vibrational  angular  fre¬ 
quency,  Uj  the  displacement  amplitudes  and  summation  con¬ 
vention  is  used  for  repeated  Cartesian  indices.  The  matrix 
M if  is  given  by 

Mn  =  cijkikjkk,  (2) 

where  ci;W  is  the  elastic  constant  tensor  elements,  and  k;  the 
wave  vector  components  of  the  vibrational  wave.  Using 
Voigt  notation,  in  the  case  of  an  hexagonal  crystal  (will  be 
explained  later),  the  Mn  matrix  reduces  to 


C\\k2x  +  C66k~  +  C44k2 

(C 12  +  Cfa)kxky 

(Cl  3  +  C44)kxkz 


(C 12  +  C66)kxky 
C66kl  +  Cnk 2  +  C44k2z 

(Ci  3  +  C44)kykz 


(Cj3  +  C44)kxkz 
(C 13  +  C44)kykz 
CUkl  +  fy  +  C^k; 


(3) 


For  any  direction  the  eigenvectors  will  be  a  quadratic 
function  of  Ai=|k|  and  hence  oj=cn(k)k  with  cn(k )  the  sound 
velocity  for  polarization  n  and  direction  k.  The  results  for 
hexagonal  crystals  were  reported  before  by  Rosen  and 
Klimker. 15  For  example,  for  k  along  [001],  kx=ky= 0  and  kz 
=  k,  the  sound  velocities  become 

crA([001])  -  \  C44/P,  with  twofold  degeneracy 

Cm([001])  =  IC^P,  (4) 

where  TA  and  LA  stand  for  transverse  and  longitudinal 
acoustic. 

Similarly,  for  any  acoustic  wave  with  wave  vector  in  the 
basal  plane,  the  sound  velocities  are  given  by  the  results  for 
the  [100]  wave  vector  direction. 


cla([100])  =  VCn/p, 


crA[oio]([iOO])  -  \  C66/p, 


c7u[ooi]([100])  -  \  C44/p,  (5) 

where  the  first  TA  mode  has  displacements  in-plane  and 
the  second  one  has  displacements  perpendicular  to  the  basal 
plane.  The  results  are  independent  of  direction  in  the  plane 
because  we  diagonalized  a  quadratic  form  representing  an 
ellipsoid.  Since  this  ellipsoid  must  have  sixfold  symmetry 
about  the  z  axis,  it  must  be  an  ellipsoid  of  revolution. 

For  a  wave  vector  of  the  form  kx=k  sin  9,  kv  =  0,  k, 
=  k  cos  9,  the  sound  velocities  are  given  by 


pc,(9)2  =  C66  sin2  6+  C44  cos2  9, 


Cn  sin2  6+  C33  cos2  6+  C44 
pc±(0)  = - — - 


(Cn  -  C44)sin2  0+  (C44-  C33) cos~  9 


+  (C13  +  C44)2cos2  9  sin2  0.  (6) 
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TABLE  I.  Calculated  equilibrium  lattice  constants  of  studied  wurtzite  semiconductors  compared  to  ex¬ 
perimental  values. 


SiC 

GaN 

InN 

ZnO 

CdSe 

a(A) 

Present 

3.05 

3.19 

3.52 

3.21 

4.28 

Expt. 

3.079“ 

3.19b 

3.54c 

3.26, d  3.253“ 

4.30, e  4.3021 

da 

Present 

1.64 

1.63 

1.61 

1.60 

1.63 

Expt. 

1.64“ 

1.627b,g 

1 ,609c 

1.60“ 

1.63e,f 

u 

Present 

0.38 

0.38 

0.38 

0.38 

0.38 

Expt. 

0.376“ 

0.382“ 

V(  A3) 

Present 

40.56 

44.59 

61.03 

45.66 

111.52 

aXRD  experiment  by  Schulz  and  Thiemann  (Ref.  22). 
bSynchrotron  EDXD  experiment  by  Xia  et  al.  (Ref.  8). 
CXRD  experiment  by  Osamura  et  al.  (Ref.  23). 
dSynchrotron  EDXD  experiment  by  Desgreniers  (Ref.  9). 
eX-ray  powder  method  by  Hotje  et  al.  (Ref.  24). 
fXRD  experiment  by  Sowa  (Ref.  25). 
gXAS  experiment  by  Perlin  et  al.  (Ref.  26). 


This  reduces  to  the  in-plane  limit  for  0=  ttI  2  and  along  c 
limit  for  0=0. 

As  will  be  seen  below,  at  zero  pressure,  we  find  for  all 
wurtzite  crystals  studied  here  that  C44<  C66.  Hence,  the  low¬ 
est  TA  sound  velocity  occurs  for  polarization  along  c  and  a 
wave  vector  in  the  plane  or  a  wave  vector  along  c  and  po¬ 
larization  in  the  plane.  Obviously  C tI  and  C33  corresponding 
to  the  longitudinal  sound  velocities  are  significantly  higher. 


energy  for  z/  =  1/2.  Here,  11  is  the  conventional  wurtzite  inter¬ 
nal  parameters.  The  structure  is  fivefold  coordinated  from 
that  point  on. 

The  b/a  distortion  can  be  written  as 


-e  0  0 

0  e  0 

0  0  0 


(8) 


B.  Relation  to  strains  and  phase  transition 

To  see  the  relation  between  specific  sound  waves  to  the 
wurtzite-to-rocksalt  phase  transitions,  we  analyze  the  strains 
related  to  the  phase  transitions  and  strains  induced  by  related 
sound  waves.  A  pure  c!a  strain  can  be  written  as 


For  e  >  0 ,  it  corresponds  to  a  compression  along  the  chosen 
b  direction  and  a  corresponding  expansion  perpendicular  to 
it  in  the  plane  and  no  distortion  along  c.  This  is  a  pure  in¬ 
plane  shear  strain  and  corresponds  to  the  limit  of  a  transverse 
sound  wave  at  45°  from  it.  In  fact,  for  a  wave  written  as 
uy(x)=u°y  exp(/k,x),  the  strain  matrix  has  the  form 


-  e/2 

0 

0 

0 

e 

0 

0 

-  e/2 

0 

(7) 

e 

0 

0 

0 

0 

e 

0 

0 

0_ 

For  e<  0,  it  corresponds  to  a  compression  along  c  combined 
with  an  expansion  in  the  plane  in  all  directions  such  that  the 
volume  is  preserved.  This  strain  maintains  the  hexagonal 
symmetry  of  the  crystal.  As  discussed,  elsewhere,  this  trans¬ 
formation  by  itself  could  transform  the  wurtzite  into  the  so- 
called  HX-structure.1'3  The  latter  differs  from  wurtzite  by  the 
fact  that  the  layers  become  unbuckled  or  fiat.  At  some  critical 
value  of  this  strain,  the  internal  coordinate  finds  a  minimum 


which  when  diagonalized,  gives  a  strain  of  the  form  of  Eq. 
(8)  with  eigenvectors  rotated  45°  from  the  original  axes  in 
the  plane.  This  means  that  the  in-plane  transverse  acoustic 
sound  waves  at  an  angle  45°  from  the  b  directions  corre¬ 
spond  exactly  to  the  symmetry  breaking  strain  involved  in 
the  wurtzite  to  rocksalt  transition,  the  strain  we  usually  call 
the  bla  strain.  Thus  the  excitation  of  such  sound  waves 
could  be  viewed  as  the  initial  step  in  creating  a  bla  strain. 


TABLE  II.  Bulk  moduls  B,  its  pressure  derivative  B'  and  elastic  constants  (in  GPa)  for  wurtzite  SiC 
(2H-SiC)  compared  with  experimental  data  on  6H-SiC. 


B 

(GPa) 

B' 

c„ 

^12 

Cl3 

C33 

C44 

Q  6 

Present 

229 

3.6 

541 

117 

61 

586 

162 

212 

Expt.“ 

220 

501 

111 

52 

553 

163 

195 

“Brillouin  scattering  experiment  by  Kamitani  et  al.  (Ref.  27). 
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TABLE  III.  Bulk  moduls  B,  its  pressure  derivative  B’  and  elastic  constants  (in  GPa)  for  wurtzite  GaN 
compared  with  other  calculations  and  experimental  data. 


B 

(GPa) 

B' 

Cu 

Cl2 

^13 

C33 

C44 

C(,6 

Present 

207 

4.2 

367 

135 

98 

409 

98 

116 

Other3 

207 

4.5 

346 

148 

105 

389 

76 

99 

Otherb 

202 

367 

135 

103 

405 

95 

116 

Other" 

201 

4.3 

366 

139 

98 

403 

97 

Expt.d 

210 

390 

145 

106 

398 

105 

122 

Expt." 

188 

3.2 

Expt.1 

245 

4.0 

Expt.g 

237 

4.3 

aDFT(LDA)  calculation  by  Kim  et  al.  (Ref.  28). 
bDFT(LDA)  calculation  by  Wright  (Ref.  16). 

CDFT(LDA)  calculation  by  Lepkowski  et  al.  (Ref.  29). 
dBrillouin  scattering  experiment  by  Polian  et  al.  (Ref.  30). 
"Synchrotron  EDXD  experiment  by  Xia  et  al.  (Ref.  8). 
fXAS  experiment  by  Perlin  et  al.  (Ref.  26). 
gXRD  experiment  by  Ueno  et  al.  (Ref.  31). 


£/=^x<Vv  (ID 

where  e ,  =  exx,  e2=eyy,  e3=ezz,  e4=2eyv  ei  =  2ezx,  and  e6 

=  2exy.  The  elastic  constants  follow  a  similar  contraction  of 
indices  rule,  e.g.,  C12 = cxxyy,  C44=Cyzyz,  etc.  Since  for  a  hex¬ 
agonal  material  there  are  5  independent  elastic  constants,  we 
can  choose  5  independent  strains,  calculate  their  elastic  en¬ 
ergy  per  unit  volume  and  write  it  in  terms  of  the  appropriate 
combination  of  elastic  constants.  We  will  obtain  5  equations 
with  5  unknowns  from  which  the  elastic  constants  can  be 
extracted.  Following  Wright  et  al.,16  we  can  use  the  follow¬ 
ing  strains  and  obtain  the  combination  of  elastic  constants 
corresponding  to  each  of  them: 

e  =  (£,£,0,0,0,0,0)->  U=(Cn  +  Cn)e1, 
e  =  (£,£,- 2f, 0,0,0) U  =  (Cu  +  C12- 4CU  + 2C33)e2, 

e  =  (0,0,  e,0,0,0)  -►  U=  -C^e2, 

In  Voigt  notation,  this  becomes  a  matrix  equation  2 


TABLE  IV.  Bulk  moduls  B,  its  pressure  derivative  B'  and  elastic  constants  (in  GPa)  for  wurtzite  InN 
compared  with  other  calculations  and  experimental  data. 


B 

(GPa) 

B' 

C„ 

Cl2 

Cl3 

C33 

C44 

Cm 

Present 

151 

4.8 

232 

115 

96 

239 

52 

59 

Other3 

147 

3.4 

220 

120 

91 

249 

36 

50 

Otherb 

141 

223 

115 

92 

224 

48 

54 

Other" 

146 

3.9 

229 

120 

95 

234 

49 

Expt.d 

125 

12.7 

3DFT(LDA)  calculation  by  Kim  et  al.  (Ref.  28). 
bDFT(LDA)  calculation  by  Wright  (Ref.  16). 

"DFT(LDA)  calculation  by  Lepkowski  et  al.  (Ref.  29). 
dXRD  experiment  by  Ueno  et  al.  (Ref.  31). 


The  elastic  constant  defining  the  sound  velocity  for  this  type 
of  acoustic  sound  waves  is  C(](v  On  the  other  hand,  sound 
waves  in  the  plane  with  displacement  along  the  c  axis  or  with 
wave  vectors  along  c  involve  C44.  Now,  because  C44<C66, 
we  can  see  that  sound  waves  producing  the  bl a  distortion  are 
not  the  easiest  to  excite.  For  any  finite  wave  vector,  they 
require  a  bit  higher  energy.  However,  as  we  will  see  at  some 
pressure  the  situation  reverses  and  the  C66  becomes  the  lower 
one.  It  is  thus  of  interest  to  check  whether  the  pressure  where 
this  crossing  happens  is  in  some  ways  related  to  the  initiation 
of  the  phase  transition. 

III.  COMPUTATIONAL  METHOD 
A.  Elastic  constant  calculation 

First,  we  briefly  recall  how  to  calculate  elastic  constants. 
The  elastic  energy  (per  unit  volume)  can  be  written  as 

^=Tey  cijkl€kl-  (10) 
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TABLE  V.  Bulk  moduls  B,  its  pressure  derivative  B'  and  elastic  constants  (in  GPa)  for  wurtzite  ZnO 
compared  with  other  calculations  and  experimental  data. 


B 

(GPa) 

B' 

Cu 

Cl2 

Cl3 

C33 

C44 

Cm 

Present 

162 

4.0 

227 

133 

118 

232 

40 

47 

Other3 

217 

117 

121 

225 

50 

50 

Otherb 

231 

111 

104 

183 

72 

60 

Expt.c 

206 

118 

211 

44.3 

44.6 

Expt.d 

181 

4.0 

209.7 

121.1 

105.1 

210.9 

42.47 

44.3 

Expt.e 

206 

4.0 

207 

117.7 

106.1 

209.5 

44.8 

44.65 

aDFT(LDA)  calculation  by  Gopal  and  Spaldin  (Ref.  32);  for  C66  we  used  the  value  from  (Cu-Ci2)/2. 
bAtomistic  simulation  techniques  based  on  the  shell  model  calculation  by  Zaoui  and  Sekkal  (Ref.  33). 
cBrillouin  scattering  experiment  Carlotti  et  al.  (Ref.  34). 

dB  and  B'  from  Synchrotron  EDXD  experiment  by  Decremps  et  al.  (Ref.  10),  C,;  from  Ultrasound  measure¬ 
ments  by  Bateman  (Ref.  35). 

eB  and  B'  from  Synchrotron  EDXD  experiment  by  Wu  et  al.  (Ref.  11),  Ctj  from  Kobiakov  (Ref.  36). 


e=  (0,0, 0,0,0,  e)  ->  U=  ^-(Cu-  Cn)<r, 


e  =  (0,0,0,  e,e,0)-^t/=C44e2.  (12) 

Note  that  for  the  calculations  under  each  of  these  strains 
the  internal  coordinates  are  allowed  to  relax.  Also,  note  that 
in  the  latter  two  cases,  the  hexagonal  symmetry  is  broken  by 
the  distortion.  Now,  we  want  to  calculate  all  these  elastic 
constants  as  functions  of  pressure.  We  first  determine  the 
undistorted  structure  total  energy  as  function  of  volume,  and 
from  p  =  -dEI  dV  determine  the  pressure  corresponding  to 
each  volume.  In  practice,  we  fit  the  Murnaghan  equation  of 
state  to  the  E{V)  results  and  invert  it  to  obtain  volume  as 
function  of  pressure. 

At  each  pressure,  and  the  corresponding  volume,  we  ap¬ 
ply  the  strains  defined  in  Eq.  (12)  and,  thus,  obtain  the  elastic 
energy  and  consequently  the  pressure  dependent  elastic  con¬ 
stants.  Although  some  of  these  strains  are  volume  conserv¬ 
ing,  e.g.,  the  second  and  last  two,  the  others  are  not  and 
contain  a  hydrostatic  component.  Still,  we  use  the  volume  at 
pressure  p  before  applying  the  distortion  in  defining  the  en¬ 
ergy  per  unit  volume  in  those  cases,  which  is  sufficient  be¬ 
cause  we  only  consider  linear  elasticity  theory.  In  other 
words,  we  make  a  quasiharmonic  approximation. 


B.  Computational  details 

To  obtain  all  total  energies  required  for  the  calculation  of 
the  elastic  constants  we  use  the  Kohn-Sham  density  func¬ 
tional  theory1718  in  the  local  density  approximation  with  the 
Perdew-Zunger  parametrization  of  exchange  and 
correlation. 19  The  total  energies  are  calculated  using  the  full- 
potential  linearized  muffin-tin  orbital  method  as  described  by 
Methfessel  and  van  Schilfgaarde.20  Well  converged  double 
basis  sets  are  utilized  and  the  Brillouin  zone  integrations  are 
done  using  a  4X4X4  Monkhorst-Pack  sampling  set.21 

IV.  RESULTS 

A.  Elastic  constants  and  sound  velocities  at  zero  pressure 

First,  we  give  our  results  obtained  for  the  equilibrium 
properties  of  the  materials  investigated  here  compared  to  ex¬ 
perimental  results.  Table  I  shows  the  wurtzite  lattice  con¬ 
stants  a,  da,  internal  parameter  u  and  the  volume  of  the  unit 
cell.  Excellent  agreement  is  obtained  in  all  cases. 

Tables  II-VI  give  our  results  for  the  elastic  constants  at 
zero  pressure  compared  with  literature  values.  Overall,  good 
agreement  is  obtained  both  with  experimental  and  calculated 
values  by  others.  For  SiC,  the  wurtzite  phase,  which  is  usu¬ 
ally  called  2H-SiC  is  difficult  to  grow  and,  therefore,  no 
experimental  data  on  elastic  constants  of  2H-SiC  are  avail- 


TABLE  VI.  Bulk  moduls  B,  its  pressure  derivative  B'  and  elastic  constants  (in  GPa)  for  wurtzite  CdSe 
compared  with  other  calculations  and  experimental  data. 


B 

(GPa) 

B' 

Cu 

Cl2 

Cl3 

C33 

C44 

C66 

Present 

60 

4.6 

80 

47 

40 

92 

15 

17 

Expt.a 

74.9 

46.09 

39.36 

84.51 

13.15 

14.40 

Expt.b 

53.3 

4.0 

74.1 

45.2 

39.3 

83.6 

13.17 

14.45 

“Ultrasound  measurements  by  Cline  et  al.  (Ref.  37). 

bB  and  B'  from  XRD  experiment  by  Sowa  (Ref.  25),  Ct]  from  resonance  measurement  by  Berlincourt  et  al. 
(Ref.  38). 
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TABLE  VII.  The  second-order  polynomials  describing  the  pres¬ 
sure  dependence  of  selected  elastic  constants  for  SiC,  GaN,  InN, 
ZnO,  and  CdSe.  The  polynomials  were  obtained  from  a  fit  to  the 
calculated  results  at  discrete  pressures.  The  pressure,  P,  is  in  GPa 
and  the  ranges  used  for  the  fit  are  shown  in  the  last  column. 


Elastic  constant 
(GPa) 

P  range 
(GPa) 

SiC 

C| !  =  542  +  4.79P-  0.0 140P2 

C33  =  586  +  8.02P-  0.0262P2 
C44=  163+0.91P-0.0015P2 
C66=213  +  0.83P-0.0038P2 

0-150 

GaN 

CU  =  369+3.74P-0.0135P2 

C33  =  405  +  4.54P-  0.0083P2 
C44=98+0.58P-0.0023P2 

C66=  1 17  +  0.48P-0.0052P2 

0-100 

InN 

CU  =  232  +  3.86P-0.017P2 

C33  =  240+ 4.72P  -  0.036P2 
C44=52  +  0.24P-0.001P2 

C66  =  58  -  0.08P- 0.005 P2 

0-40 

ZnO 

Cn  =  227  +  1 .62P  +  0.073P2 
C33=232+  1.16P+0.165P2 
C44=40+0.03P-0.005P2 
C66=47-0.86P +  0.006P2 

0-16 

CdSe 

C„  =  81  +  3.53P-0.075P2 

C33  =  91+5.49P-0.197P2 

C44=  15+0.05P-0.026P2 
C66=17  +  0.11P-0.111P2 

0-10 

able.  However,  one  finds  that  the  elastic  constants  for  4H  and 
6H  SiC,  two  slightly  different  polytypes,  are  very  close  to 
each  other.  We  thus  compare  with  the  values  for  6H. 

The  sound  velocities  at  zero  pressure  extracted  from  these 
at  zero  pressure  are  given  in  Tables  VIII  for  wave  vector 
directions  along  [001]  and  [100].  As  discussed  above,  the 
sound  velocities  for  wave  vectors  in  the  plane  are  indepen¬ 
dent  of  the  direction  angle  in  the  plane.  For  wave  vectors  that 
are  intermediate  between  in-plane  and  perpendicular  to  the 
basal  plane,  the  sound  velocities  have  a  well  defined  angular 
dependence  given  in  Eq.  (6). 

B.  Pressure  dependent  elastic  constants  and  sound  velocities 

In  Figs.  1  and  2,  we  show  the  pressure  dependence  of  the 
elastic  constants  and  sound  velocities,  respectively.  The  dots 
show  the  first  principles  results  for  the  given  pressures.  The 
lines  are  the  second-order  polynomial  fit  to  the  results.  The 
fitting  results  (quadratic  equations)  for  the  elastic  constants 
are  shown  in  Table  VII.  We  can  see  quite  different  behaviors 
depending  on  the  material  and  the  elastic  constants  consid¬ 
ered.  The  compressional  moduli,  Cn  and  C33  are  increasing 
with  pressure  but  for  most  materials  we  see  that  C33  curves 
downward  at  higher  pressures.  The  exception  is  ZnO,  where 
it  curves  upward.  For  SiC  and  GaN,  the  shear  type  elastic 
constant  C66  initially  increases  but  goes  through  a  maximum 
and  then  bends  down.  For  InN,  ZnO,  and  CdSe,  it  monotoni- 


cally  decreases  with  increasing  pressure  starting  from  zero 
pressure.  In  ZnO,  it  slightly  curves  upward  instead  of  down¬ 
ward.  The  C 44  changes  from  increasing  with  pressure  in  SiC, 
GaN,  and  InN,  to  staying  almost  constant  in  ZnO  and  de¬ 
creasing  in  CdSe.  Since  the  sound  velocities  are  directly  de¬ 
rived  from  these  elastic  constants,  similar  behavior  is  also 
observed.  Their  behavior  is,  however,  slightly  different  be¬ 
cause  the  mass  density  involved  in  the  sound  velocity  also 
changes  as  a  function  of  pressure.  To  the  best  of  our  knowl¬ 
edge,  no  measurements  of  the  sound  velocities  or  elastic  con¬ 
stants  as  function  of  pressure  are  available  for  these  materi¬ 
als.  Our  calculations  predict  and  interesting  variety  of 
behaviors,  which  would  be  beneficial  to  verify  experimen¬ 
tally.  The  elastic  constants  as  function  of  pressure  were  pre¬ 
viously  calculated  using  DFT  by  tepkowski  et  air’  for  GaN 
and  InN.  They  extracted  the  linear  and  quadratic  pressure 
coefficients.  Their  results  for  the  two  materials  are  compa¬ 
rable  to  ours  (as  shown  in  Fig.  1),  but  we  considered  a  larger 
pressure  range.  For  ZnO,  the  elastic  constants  as  function  of 
pressure  were  previously  calculated  by  Zaoui  and  Sekkal33 
but  using  only  pair  potential  shell  model  level  of  calcula¬ 
tions.  For  CdSe,  the  elastic  constants  as  a  function  of  pres¬ 
sure  were  previously  studied  for  only  the  zincblende  phase 
by  Deligoz  et  al .41 

C.  Relation  to  phase  transition 

As  discussed  in  the  introduction,  it  is  of  interest  to  con¬ 
sider  the  pressure  where  C44  crosses  with  C66.  Above  this 
pressure,  transverse  sound  waves  in  the  basal  plane  should 
become  easier  to  excite  than  sound  waves  related  to  da 
strains.  The  transverse  sound  waves  are  closely  related  to  the 
symmetry  breaking  bl  a  strain  leading  to  the  phase  transition, 
while  the  sound  waves  related  to  da  strains  preserve  the 
hexagonal  symmetry.  In  order  to  initiate  the  phase  transition, 
it  might  seem  important  to  break  the  symmetry. 

We  thus  consider  the  pressures  where  these  two  shear 
elastic  constants  cross,  or  equivalently  where  the  sound  ve¬ 
locities  cross.  The  values  are  summarized  in  Table  IX  along 
with  the  equilibrium  phase  transition  pressure  (the  pressure 
where  the  calculated  enthalpies  of  the  rocksalt  and  wurtzite 
phases  are  equal)  and  the  experimentally  observed  phase 
transition  pressures.  We  estimate  the  uncertainty  on  these 
crossover  pressures  to  be  of  the  order  a  few  GPa. 

We  can  see  that  for  SiC,  GaN,  and  CdSe  the  cross-over 
pressure  is  higher  than  the  actual  transition  pressure,  while 
the  equilibrium  transition  pressure  is  lower  than  the  experi¬ 
mental  value.  The  latter  has  been  previously  attributed  to  the 
existence  of  a  barrier  between  the  two  phases  at  the  transition 
pressure.  It  can  be  considered  a  kinetic  effect.  For  InN,  our 
calculated  equilibrium  transition  pressure  is  in  good  agree¬ 
ment  with  the  experimental  value  while  the  crossover  pres¬ 
sure  is  higher.  Note  that,  generally  the  InN  crystal  quality  is 
quite  poor  so  the  measurement  of  this  particular  material 
might  contain  a  sizable  error  bar.  For  ZnO,  the  calculated 
crossover  pressure  is  slightly  higher  than  the  calculated  equi¬ 
librium  transition  pressure.  The  observed  experimental  tran¬ 
sition  pressures  of  9.1-10.5  GPa  exceeds  both  of  these  val¬ 
ues  but  by  only  ~1  GPa,  which  is  well  within  the 
computational  error  bar. 
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TABLE  VIII.  Sound  velocities  (in  km/s)  at  zero  pressure  for  SiC,  GaN,  InN,  ZnO,  and  CdSe  calculated 
from  the  elastic  constants,  k  is  the  wave  vector  and  n  is  the  polarization.  The  Phillips  ionicity  parameters  (f ,) 
are  also  listed. 


fi 

k-> 

n— > 

[100] 

[001] 

[100] 

[001] 

[010] 

[001] 

[100] 

SiC 

0.177 

Present 

12.8 

7.0 

8.0 

13.4 

7.0 

Other* 

12.5 

7.1 

7.8 

13.1 

7.1 

Expt.b 

12.21 

7.69 

Expt.c 

13.27 

7.24 

GaN 

0.500 

Present 

7.67 

3.96 

4.32 

8.10 

3.96 

Other4 

8.0 

4.1 

6.3 

8.0 

4.1 

InN 

0.578 

Present 

5.76 

2.73 

2.89 

5.84 

2.73 

Other4 

5.3 

1.2 

2.5 

5.2 

1.2 

ZnO 

0.616 

Present 

6.23 

2.61 

2.83 

6.30 

2.61 

Other** 

6.08 

2.92 

2.92 

6.19 

2.92 

Expt.f 

6.08 

2.73 

2.79 

6.10 

2.73 

CdSe 

0.699 

Present 

3.76 

1.62 

1.71 

4.02 

1.62 

Expt.g 

3.63 

1.52 

1.59 

3.86 

1.52 

“Calculated  from  elastic  constants  by  Kamitani  et  al.  (Ref.  27). 
bUltrasound  measurements  by  Schreiber  et  al.  (Ref.  39). 
cRaman  measurements  on  various  polytypes,  Feldman  et  al.  (Ref.  40). 
dCalculated  from  elastic  constants  by  Wright  et  al.  (Ref.  16). 
“Calculated  from  elastic  constants  by  Gopal  and  Spaldin  (Ref.  32). 
TJltrasound  measurements  by  Bateman  (Ref.  35). 
gUltrasound  measurements  by  Cline  et  al.  (Ref.  37). 


As  discussed  before  by  Limpijumnong  and  Jungthawan,42 
for  more  ionic  or  softer  tetrahedral  semiconductors,  the  cal¬ 
culated  transition  pressures  are  close  to  the  experimental 
ones,  while  for  the  stiffer  ones,  the  calculated  equilibrium 
transition  pressures  are  significantly  underestimated.  The  un¬ 
derestimation  of  the  calculated  transition  pressure  was  attrib¬ 
uted  to  the  kinetic  barrier  in  the  actual  transition.  Here,  we 
observe  some  relationship  between  the  crossover  pressure 
and  the  transition  pressure.  We  find  that  the  crossover  pres¬ 
sure  might  serve  as  an  upper  limit  of  the  transition  pressure. 
In  addition,  we  propose  that  in  combination  with  the  calcu¬ 
lated  equilibrium  transition  pressure  it  can  be  used  to  im¬ 
prove  the  prediction  of  the  actual  transition  pressure,  as  will 
be  described  later.  This  observation  is  especially  useful  for 
the  stiffer  materials,  for  e.g.,  SiC  and  GaN  in  our  present 
study,  where  the  calculated  transition  pressure  is  significantly 
underestimated. 

Above  the  crossover  pressure,  the  transverse  sound  waves 
associated  with  C66  turn  easier  to  excite  than  sound  waves 
associated  with  C44.  However,  even  before  the  C66  and  C44 
cross,  one  may  expect  that  as  they  approach  each  other  it 
already  become  probable  to  excite  the  transverse  acoustic 
waves  at  finite  temperature  and  ultimately  initiate  the  transi¬ 
tion.  Note  that  in  this  model,  we  consider  the  excitation  of 
such  modes  as  the  triggering  effect.  Once  the  system  starts 
exhibiting  a  strain  of  this  type,  the  da  distortion  follows  it 
through  the  anharmonic  interactions  between  these  two  types 
of  modes.  This  then  leads  the  system  toward  the  barrier  tran¬ 
sition  point  which  separates  the  wurtzite  from  the  rocksalt 


valley  as  can  be  seen  from  the  energy  landscapes  given  in 
Sarasamak  et  al .3  Note  that  the  excitation  energy  of  the  pho¬ 
non  with  the  wavelength  of  about  1  nm  can  be  estimated 
from  hw=hvk~  10-20  meV,  which  is  comparable  to  the 
thermal  energy  at  room  temperature.  This  means  that  within 
a  nucleation  region  of  about  1  nm  the  homogeneous  trans¬ 
formation  might  be  triggered  already  at  room  temperature 
and  serves  as  a  nucleus  of  the  new  phase. 

From  Table  IX,  it  is  clear  that  the  calculated  equilibrium 
transition  pressure  p,  generally  underestimates  the  experi¬ 
mentally  observed  pressure  pu  by  about  as  much  percentage 
wise  as  the  crossover  pressure  overestimates  it.  Therefore, 
we  propose  that  a  good  and  simple  approximation  to  the 
experimental  pressure  would  be  ( p,+pc)/2 .  Our  calculated 
values  of  (/?,+/?,)/ 2  are  tabulated  in  Table  IX  and  are  very 
close  the  experimental  values  for  SiC  and  GaN.  As  a  further 
test  that  this  gives  a  useful  estimate  in  particular  for  cases 
where  there  is  evidence  of  a  significant  transition  kinetic  bar¬ 
rier,  we  consider  AIN.  Based,  on  the  data  of  Lepkowski  et 
al.29  for  the  pressure-dependent  elastic  constants,  we  obtain 
the  crossover  pressure  to  be  39.6  GPa.  The  calculated  equi¬ 
librium  transition  pressure  by  Serrano  et  a/.43  is  9.2  GPa. 
This  gives  as  estimate  for  the  actual  transition  pressure  {pt 
+/?,.) /2  =  24.4  GPa.  The  experimental  value  by  Ueno  et  al.44 
is  22.9  GPa. 

Note  that  the  correlation  found  here  between  the  sound 
wave  excitations  of  a  particular  symmetry  and  the  transition 
pressures  is  suggestive  but  not  a  complete  proof  of  their  roles 
in  the  actual  transition  processes.  The  actual  transition  re- 
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FIG.  1.  (Color  online)  Calculated  pressure  dependence  of  the 
elastic  constants  in  (a)  SiC,  (b)  GaN,  (c)  InN,  (d)  ZnO,  and  (e) 
CdSe.  In  (b)  and  (c),  the  square  symbols  show  the  calculated  results 
by  Lepkowski  et  al.  (Ref.  29). 

quires  the  occurrence  of  both  in-plane  and  perpendicular  to 
the  plane  types  of  strain  to  reach  the  transition  barrier.  The 
transition  barrier  height42  and  the  overall  enthalpy  landscape3 
therefore  play  a  significant  role  in  the  phase  transformation 
process. 

D.  Trends  with  ionicity 

Because  we  have  systematically  calculated  several  wurtz- 
ite  crystals,  i.e.,  SiC,  GaN,  InN,  ZnO,  and  CdSe,  that 
spanned  a  wide  range  of  ionicity,  we  can  study  the  trend  of 
the  elastic  constants  and  sound  velocities  with  the  material’s 
ionicity.  To  quantify  the  ionicity  of  these  materials,  we 
choose  to  use  Phillips’  ionic  scale45  (/))  which  has  the  range 
between  0  (completely  covalent)  and  1  (completely  ionic). 
This  scale  is  based  on  optical  properties  of  the  semiconduc¬ 
tors  and  has  been  shown  in  the  past  to  be  related  to  the 
preference  for  octahedral  and  tetrahedral  bonding  and  the 
transition  pressures.46-47  The  values  of  /,■  for  the  materials 
studied  are  listed  in  Table  VIII.  In  Fig.  3(a),  we  plotted  the 
elastic  constants  as  a  function  of  /,  at  zero  pressure.  We  can 
see  that  all  four  elastic  constants  monotonically  decrease  as 
fi  increases.  The  same  holds  true  at  finite  pressures.  For  the 
sound  velocity,  the  general  trend  also  follows  that  of  the 
elastic  constant  because  the  sound  velocity  is  proportional  to 
the  square  root  of  the  corresponding  elastic  constant  as 


FIG.  2.  (Color  online)  Calculated  pressure  dependence  of  the 
sound  velocities  in  (a)  SiC,  (b)  GaN,  (c)  InN,  (d)  ZnO,  and  (e) 
CdSe. 

shown  in  Eq.  (4)  and  (5).  In  Fig.  3(b),  the  sound  velocities 
are  shown  as  a  function  of  /,-.  We  note  that  in  both  the  LA 
sound  velocities  and  compressional  elastic  constants,  InN 
dips  a  bit  below  the  trend  lines.  This  may  be  indicative  of  the 

TABLE  IX.  Crossover  pressure  pc  where  C44=C66,  calculated 
equilibrium  phase  transition  pressure  p,  (taken  from  Ref.  3,  except 
for  AIN)  and  experimental  transition  pressures  pu  for  wurtzite-to- 
rocksalt  transition,  all  in  GPa,  for  SiC,  GaN,  AIN,  InN,  ZnO,  and 
CdSe. 


Pc 

Pt 

(p,+Pc)/2 

Pi, 

SiC 

131.1 

64.9 

98.0 

© 

o 

GaN 

65.5 

44.1 

54.8 

52. 2b 

AIN 

39.6C 

9.2d 

24.4 

22. 9e 

InN 

15.7 

12.2 

13.9 

12. lb 

ZnO 

8.8 

8.2 

8.4 

9.1f,  10. 5g 

CdSe 

5.2 

2.2 

3.7 

4.0h 

aXRD  experiment  by  Yoshida  et  al.  (Ref.  50). 
bXRD  experiment  by  Ueno  et  al.  (Ref.  31). 
cCalculated  here  from  data  of  Lepkowski  et  al.  (Ref.  29). 
dFrom  Serrano  et  al.  (Ref.  43). 
eFrom  Ueno  et  al.  (Ref.  44). 

Synchrotron  XRD  experiment  by  Desgreniers  (Ref.  9). 
BSynchrotron  EDXD  experiment  by  Kumar  et  al.  (Ref.  51). 
hSynchrotron  XRD  experiment  by  Wang  et  al.  (Ref.  52). 
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FIG.  3.  (Color  online)  Calculated  (a)  elastic  constants  and  (b) 
sound  velocities  as  a  function  of  ionicity  (Phillips’  scale)  at  zero 
pressure  (P= 0  GPa). 

fact  that  the  ionicity  of  InN  is  underestimated  by  the  Phillips 
scale.  InN  has  a  larger  ionicity  in  the  Garcia  scale.48  The 
ionicities  of  the  Ill-nitrides  were  discussed  in  Ref.  49 

V.  CONCLUSIONS 

We  presented  the  calculated  elastic  constants  and  sound 
velocities  in  wurtzite  SiC,  GaN,  InN,  ZnO,  and  CdSe  as 
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functions  of  pressure.  We  found  interesting  nonlinear  behav¬ 
ior  of  the  elastic  constants  and  sound  velocities  with  pressure 
that  is  quite  different  in  the  different  materials.  We  also 
showed  that  the  elastic  constants  are  lower  in  materials  with 
higher  ionicity.  Generally,  the  sound  velocities  follow  the 
same  trend  as  the  elastic  constants,  i.e.,  lower  in  materials 
with  higher  ionicity.  The  pressure  at  which  the  C44  and  C66 
shear  elastic  constants  cross  are  investigated.  When  C66  be¬ 
comes  lower,  it  should  be  easier  to  excite  the  symmetry 
breaking  strain  mode  related  to  the  wurtzite-to-rocksalt  phase 
transition.  We  found  that  the  crossover  pressure  is  generally 
higher  than  the  experimental  transition  pressure  and  ex¬ 
plained  it  by  the  fact  that  finite  temperature  fluctuations 
make  it  possible  to  excite  the  required  modes  even  before  the 
two  sound  velocities  cross.  We  proposed  the  average  of  the 
crossover  and  calculated  equilibrium  transition  pressures  to 
be  a  good  estimate  for  the  experimental  transition  pressure  in 
cases  where  there  is  a  significant  transition  enthalpy  barrier. 
The  results  suggest  that  the  in-plane  transverse  sound  waves 
could  play  a  significant  role  in  triggering  the  phase  transi¬ 
tion.  Of  course,  the  details  near  the  transition  point,  the 
height  of  the  barrier,  and  the  coupling  of  the  in  and  out  of 
plane  acoustic  modes,  which  are  both  required  to  reach  the 
transition  point,  are  also  important  for  the  transition  to  take 
place. 
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To  investigate  the  local  structure  of  InN,  ln203  and  their  alloys,  synchrotron  (In  L3-edge)  X-ray  absorption 
near  edge  structures  (XANES)  of  indium  oxynitride  samples  with  varied  0  contents  are  used  in  conjunc¬ 
tion  with  first-principles  calculations.  A  good  agreement  between  the  measured  and  simulated  spectra  is 
obtained.  It  is  found  that  the  spectra  are  sensitive  to  the  coordination  number  of  the  In  atoms,  i.e.,  four¬ 
fold  for  InN-like  structures  and  sixfold  for  ln203-like  structures.  Moreover,  the  spectra  are  quite  insensi¬ 
tive  to  the  species  (N  or  0)  around  In.  The  calculated  band  structures  and  density  of  states  are  also 
presented  and  discussed. 

©  2010  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

Indium  oxynitride  is  an  alloy  between  indium  nitride  (InN)  and 
indium  oxide  (In203).  The  crystal  structures  of  InN  and  ln203  are 
very  different,  i.e.,  wurtzite  for  InN  and  bixbyite  for  ln203.  The  band 
gap  of  InN  has  been  mistaken  to  be  ~1.9  eV  for  a  long  time  (see  the 
discussion  in  Refs.  [1,2]).  Recently,  a  high  quality  InN  has  been  ob¬ 
tained  by  molecular  beam  epitaxy  (MBE)  and  the  actual  InN  band- 
gap  of  ~0.7  eV  has  been  observed  [3,4],  In  addition,  the  bandgap  of 
ln203  remains  contentious  with  the  reported  values  ranging  from 
~3.6eV  [5]  to  below  3.0  eV  [6].  The  nitrogen  incorporation  in 
ln203  can  reduce  the  bandgap  into  the  visible  range,  which  will 
be  beneficial  for  solar  photocatalytic  applications  such  as  hydrogen 
production  from  water  [7].  Because  bandgap  of  indium  oxynitride 
can  vary  in  a  very  wide  range,  this  class  of  alloys  is  widely  used  for 
optical  coating  applications  despite  its  poor  crystal  quality,  espe¬ 
cially  when  it  is  grown  by  traditional  techniques  such  as  RF  mag¬ 
netron  sputtering  [8]. 

In  alloys  with  low  oxygen  content,  0  atoms  can  substitute  on 
the  N  site  in  the  wurtzite  InN  crystal  structure  [9].  For  alloys  with 
high  oxygen  content,  the  crystal  structure  is  still  unclear.  Note  that 
the  O  content  that  can  substitute  on  the  N  site  in  the  wurtzite 
structure  would  depend  strongly  on  the  growth  techniques  and 
conditions.  In  the  past,  N  K-edge  XANES  measurement  has  been 
performed  to  study  the  crystal  structure  of  this  class  of  alloys. 
Unfortunately,  probing  the  local  structure  of  anions  in  this  alloy 
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system  is  not  very  useful  because  the  local  structure  of  the  anion 
in  both  wurtzite  and  bixbyite  structures  are  fourfold  coordinated. 
On  the  other  hand,  the  local  geometry  of  In  in  wurtzite  and  bixby¬ 
ite  structures  are  different,  i.e.,  fourfold  in  InN  and  sixfold  in  ln203. 
As  we  have  previously  illustrated  [10],  probing  the  local  structure 
of  In  in  this  alloy  system  can  give  useful  information  on  the  local 
structure  of  the  alloys.  In  our  published  Letter  [10],  we  employed 
In  L3-edge  XANES  spectra  in  conjunction  with  ab  initio  XANES  sim¬ 
ulation  to  characterize  the  indium  oxynitride  thin  films  with  varied 
compositions.  We  have  shown  that  XANES  is  a  powerful  tool  to 
determine  the  local  structure  of  indium  oxynitride  films.  Here, 
the  detailed  results  that  have  not  been  presented  in  our  Letter 
due  to  limited  space  and  additional  results  on  the  simulation  of 
In  surrounded  by  both  oxygen  and  nitrogen  neighbors  (mixed 
neighbors)  are  presented. 


2.  Methods 

2.3.  Experimental  details 

Indium  oxynitride  films  were  grown  by  RF  magnetron  sputter¬ 
ing  at  room  temperature  using  a  technique  called  reactive  gas  tim¬ 
ing  [11,12].  Before  the  growth  process,  the  chamber  was  pre¬ 
evacuated  to  ~1CT5  Pa.  After  that,  the  N2  and  02  gas  were  flown 
interchangeably  at  the  flow  rate  of  10  standard  cm3  per  minute 
(seem)  onto  the  99.999%  purity  In  target.  The  sputtering  gas  pres¬ 
sures  were  set  at  0.34  Pa  and  0.32  Pa  for  N2  and  02,  respectively. 
The  indium  oxynitride  samples  with  varied  O  and  N  contents  can 
be  obtained  by  controlling  N2  and  02  gas  timing  [10].  The  optical 
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bandgap  and  the  0:ln  composition  of  all  samples  were  measured 
by  UV-visible  and  Auger  electron  spectroscopy  (AES),  respectively 
[10]. 

The  In  L3-edge  XANES  measurements  [10]  were  performed  in 
the  fluorescent  mode  with  a  13-component  Ge  detector  (Canber- 
a)  at  the  X-ray  absorption  spectroscopy  beamline  (BL-8)  of  the 
Siam  Photon  Source  (electron  energy  of  1.2  GeV,  beam  current 
120-80  mA),  Synchrotron  Light  Research  Institute,  Thailand.  A 
double  crystal  monochromator  with  Si  (1  1  1)  was  used  to  scan 
the  synchrotron  X-ray  with  the  photon  energy  step  of  0.25  eV 
in  the  range  of  3640-3950  eV,  covering  the  XANES  region  of  In 
L3-edge. 


2.2.  Computational  details 

In  order  to  understand  the  local  microscopic  structure  of  in¬ 
dium  oxynitride,  first-principles  XANES  simulations  of  wurtzite 
InN,  bixbyite  ln203,  and  different  alloy  models  of  (InN)x(InO0.5)i_x 
were  performed.  The  detailed  crystal  geometries  were  optimized 
based  on  first-principles  density  functional  calculations.  The  den¬ 
sity  functional  theory  with  local  density  approximation  (LDA) 
and  ultrasoft  Vanderbilt  pseudo  potentials  (USPP)  [13]  as  imple¬ 
mented  in  the  VASP  codes  were  used  [14,15].  The  cutoff  energy 
for  the  plane  wave  basis  set  was  400  eV.  For  fc-point  sampling, 
Monkhorst-Pack  scheme  [16]  was  used  for  Brillouin  zone  integra- 


Fig.  1.  The  local  structures  around  In  atoms  in  (a)  pure  bixbyite  ln203,  (b)  pure  wurtzite  InN,  (c)  wurtzite  InN  with  four  surrounding  N  atoms  replaced  by  O  atoms,  and  (d) 
wurtzite  InN  with  one  surrounding  N  atom  replaced  by  an  O  atom  (left)  and  bixbyite  ln203  with  one  O  atom  replaced  by  a  N  atom  (right).  All  bond  distances  are  given  as  a 
percentage  difference  from  an  average  ln203  bond  distance  (dcaic  =  2.170  A). 
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tion  (7  x  7  x  7  for  four-atom  wurtzite  cell,  2  x  2  x  2  for  72-atom 
wurtzite  supercell,  and  3x3x3  for  40-atom  bixbyite  cell).  All 
atoms  were  allowed  to  relax  until  the  residue  forces  acting  on  each 
atom  were  less  than  10~3  eV/A.  The  volumes  of  the  cells  used  in  the 
calculations  are  based  on  the  fully  relaxed  lattice  constant  of  bulk 
wurtzite  InN  (or  bixbyite  ln203).  The  relaxed  local  structures  sur¬ 
rounding  In  atoms  in  InN,  ln203,  and  selected  InNx03_x  alloy  mod¬ 
els  are  shown  in  Fig.  1.  All  In  atoms  in  wurtzite  InN  crystal  are 
equivalent.  However,  there  are  two  species  of  In  atoms  in  ln203  (la¬ 
beled  Ini  and  In2)  with  the  composition  ratio  of  1 :3.  For  the  alloys, 
a  few  selected  local  structure  models  surrounding  an  In  atom  were 
tested.  For  the  first  model,  we  replaced  four  N  atoms  surrounding 
an  In  atom  with  0  atoms  in  a  wurtzite-InN  supercell.  For  the  sec¬ 
ond  model,  only  one  N  atom  was  replaced  by  an  0  atom  in  a  wurtz- 
ite-InN  supercell.  For  the  third  model,  only  one  0  atom  was 
replaced  by  a  N  atom  in  a  bixbyite  ln203  cell.  These  tests  were  per¬ 
formed  to  see  the  effect  of  neighboring  specie  on  the  simulated 
XANES  spectrum.  Note  that,  to  obtain  the  proper  local  geometries, 
after  the  replacement  of  N  by  0  (or  0  by  N  in  the  third  model),  all 
atoms  in  the  cell  were  allowed  to  relax  by  first-principles  calcula¬ 
tions  before  they  were  further  used  in  XANES  calculations. 

To  simulate  ab  initio  XANES  of  InN,  ln203,  and  InNx03_x  alloys, 
we  used  the  FEFF8.2  codes  [17,18].  The  codes  utilize  a  full  multiple 
scattering  approach  based  on  ab  initio  overlapping  muffin-tin 
potentials.  The  muffin-tin  potentials  used  in  FEFF  codes  are  ob¬ 
tained  using  self-consistent  calculations  with  Hedin-Lundqvist  ex- 
change-correlation  function  [19].  The  self-consistent  calculations 
were  performed  in  the  sphere  radius  4  A  (containing  approxi¬ 
mately  40-atoms)  around  the  absorber  In  atom.  The  full  multiple 
scattering  calculations  include  all  possible  paths  within  a  larger 
cluster  radius  of  7.4  A  (containing  approximately  140  atoms)  [10]. 

The  electronic  transitions  associated  with  the  XANES  measure¬ 
ments  must  follow  the  dipole  selection  rule.  Consequently,  i3-edge 
XANES  corresponds  to  the  (s  +  d)-partial  density  of  state  (PDOS).  As 
a  result,  PDOS  can  also  be  used  to  compare  with  the  XANES  mea¬ 


surement  [20].  An  X-ray  absorbance  p(co)  is  given  by  the  Fermi’s 
golden  rule, 

Moc^Kfim^iEi-Ef  +  co),  (1) 

f 

where  |i),  |/),  £,-,  and  £/are  the  initial  and  final  states  and  their  ener¬ 
gies,  respectively,  to  and  D  are  the  photon  frequency  and  dipole 
operator. 

3.  Results  and  discussion 

3.3.  Electronic  structures  of  InN,  ln203,  and  InNi_xOx  alloy 

The  calculated  cell  parameters  of  wurtzite  InN  are  a  =  3.51  A 
and  c/a  =  1.60  A  (these  are  to  compare  with  the  experimental  val¬ 
ues  of  3.533  A  and  1.611  A,  respectively  [21]).  We  obtained  a 
slightly  negative  bandgap  using  first-principles  calculations  as 
illustrated  in  Fig.  2a.  The  density  of  states  (DOS)  are  shown  in 
the  top  panel  of  Fig.  3  with  site  decomposed  partial  DOS  of  In 
and  N  in  the  middle  panel  and  bottom  panel,  respectively.  From 
PDOS,  it  is  clear  that  the  lowest  narrow  band  centered  at  -13  eV 
is  mainly  composed  of  N  s  and  the  higher  valence  band  (VB)  states 
located  at  around  -5  to  0  eV  are  composed  of  In  s,  p  and  N  p  with  a 
small  degree  of  In  d.  The  lower  conduction  band  (CB)  states  located 
at  0-12  eV  are  the  mixing  of  In  p  and  N  p  with  some  degree  of  In  s 
and  N  s,  whereas  the  higher  CB  states  (above  12  eV)  are  mainly  In  d. 

For  ln203,  the  bixbyite  unit  cell  is  composed  of  40-atom  with 
the  calculated  lattice  constant  of  ~10.07  A  (this  is  to  compare  with 
the  experimental  value  of  10.117  A  [22]).  We  obtained  the  calcu¬ 
lated  bandgap  of  1.26  eV  which  is  in  agreement  with  the  calcula¬ 
tion  value  of  1.21  eV  obtained  by  Erhart  [23].  The  calculated 
bandgaps  of  InN  and  ln203  are  lower  than  those  of  experimental 
values  because  of  the  well  know  DFT  problems.  The  band  structure 
and  the  DOS  are  shown  in  Figs.  2b  and  4,  respectively.  The  symme¬ 
try  points  labeled  in  the  band  structure  plot  is  used  according  to 


Fig.  2.  Band  structures  of  (a)  wurtzite  InN  and  (b)  bixbyite  ln203.  The  energy  is  referenced  to  the  Fermi  level. 
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Fig.  3.  The  calculated  density  of  states  (DOS)  of  wurtzite  InN:  (top  panel)  the  total 
DOS,  (middle  and  bottom  panels)  the  site  and  angular  momentum  projected  partial 
density  of  states  (PDOS)  of  In  and  N  atoms. 


Fig.  4.  The  calculated  DOS  of  bixbyite  ln203:  (top  panel)  the  total  DOS,  (middle  and 
bottom  panels)  the  site  and  angular  momentum  projected  PDOS  of  In  and  0  atoms. 


Ref.  [23].  Similar  to  the  case  of  InN,  the  lowest  narrow  band  cen¬ 
tered  at  -17  eV  is  mainly  0  s  and  the  higher  VB  located  around 
-6  to  0  eV  are  mainly  composed  of  0  p  with  a  small  degree  of  In 
s,  p  and  d.  The  lowest  CB  states  located  at  0-7  eV  are  mainly  com¬ 
posed  of  0  p  hybridized  with  In  s  and  p  with  a  small  degree  of  In  d. 
The  higher  CB  states  (7-14  eV)  are  the  mixing  of  Op  and  In  p  with  a 
small  degree  of  In  d.  The  higher  states  (above  14  eV)  are  mainly  In 
d. 

For  the  alloy,  we  calculated  the  DOS  of  wurtzite  InN^O*  alloy 
to  determine  the  effects  of  0  incorporation.  Starting  with  a  72- 
atom  wurtzite  supercell  of  InN,  we  randomly  replaced  14  (out  of 
36)  nitrogen  atoms  with  oxygen  atoms  without  any  site  discrimi¬ 
nation.  This  gives  approximately  40%  0  content  and  the  alloy  will 
be  referred  to  as  InN0.600.4.  In  addition,  we  also  calculated  20%  oxy¬ 
gen  content  by  replaced  7  nitrogen  atoms  with  oxygen  atoms.  All 
atoms  in  the  supercell  were  allowed  to  relax  before  calculating 
the  DOS.  Because  a  considerable  fraction  of  N  atoms  were  replaced 
by  0  atoms,  the  CB  is  partially  occupied.  To  properly  simulate 
XANES,  only  states  above  the  Fermi  level  can  contribute  as  a  final 
state  of  the  excitations.  The  DOS  along  with  PDOS  are  shown  in 


Fig.  5.  The  calculated  DOS  of  InN0.600.4  alloy:  (top  panel)  the  total  DOS,  (bottom 
three  panels)  the  site  and  angular  momentum  projected  PDOS  of  In,  N,  and  0  atoms. 


Fig.  5.  It  is  clear  that  the  DOS  of  the  alloy  is  a  mixture  of  InN  DOS 
and  ln203  DOS.  The  two  lowest  narrow  bands  centered  at  -17  eV 
and  -13  eV  are  0  s  and  N  s,  respectively.  The  higher  VB  states  lo¬ 
cated  at  around  -6  to  0  eV  are  the  mixing  of  N  p,  0  p,  In  s  and  In 
p.  Again,  similar  to  the  case  of  InN  and  ln203,  the  lowest  CB  states 
located  at  0-6  eV  are  mainly  composed  of  0  p  and  N  p  hybridized 
with  In  s  and  p  with  a  small  degree  of  In  d.  The  higher  CB  states  (6- 
14  eV)  are  the  hybridization  of  N  p  (or  0  p)  with  In  (mainly  p  and  d) 
with  a  small  degree  of  In  s  and  N  s  (or  0  s).  The  higher  states  (above 
14  eV)  are  mainly  In  d. 

To  compare  the  electronic  states  obtained  from  DFT  method 
with  In  i3-edge  XANES  spectra,  we  calculated  In  s  +  d  PDOS  and  ap¬ 
plied  some  broadening.  The  simulated  spectra  of  InN  and  InNo.600.4 
are  shown  in  Fig.  6.  It  is  not  surprising  that  both  spectra  appear  to 
be  very  similar.  This  is  because  the  In  atoms  of  both  spectra  are 
fourfolded.  The  main  difference  between  the  two  simulated  spectra 
is  the  reduction  in  the  magnitude  of  the  first  peak  located  at 
around  10  eV  (marked  by  an  arrow  in  Fig.  6)  above  the  Fermi  en¬ 
ergy.  This  reduction  is  expected  because  parts  of  the  CB  of 
InN0.600.4  alloy  are  occupied  and  cannot  serve  as  final  states  of 
the  absorption.  Based  on  these  PDOS,  we  expected  that  the  i3-edge 


Fig.  6.  The  angular  momentum  (s  +  d)  projected  PDOS  of  In  atoms  of  InN  and 
lnN0600.4  alloy  in  the  wurtzite  structures.  The  dotted  lines  show  the  PDOS.  The  full 
thick  lines  show  the  smeared  CB  (s  +  d)  states. 
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spectrum  of  wurtzite-InN06Oo.4  remains  similar  to  that  of  wurtzite 
InN  except  the  reduction  in  the  first  peak.  This  reduction  of  the  first 
peak  is  in  a  good  agreement  with  the  measured  and  simulated 
XANES  spectra  (based  on  the  FEFF  codes)  that  we  previously  pre¬ 
sented  [10]  and  will  be  briefly  discussed  again  in  the  next  section. 

3.2.  XANES  spectra  of  InN,  ln203,  and  InN,.xOx  alloy 

In  our  previous  Letter  [10],  the  In  i3-edge  XANES  were  used  to 
study  the  local  structures  of  indium  oxynitride  thin  films  grown  un¬ 
der  different  conditions.  For  completeness,  the  main  content  of  the 
previous  work  will  be  briefly  given  here.  Five  samples  of  indium 
oxynitride  alloys  (N30/00,  N30/05,  N30/010,  N30/020  and  NO/ 
030,  arranged  from  low  to  high  0  content)  have  been  characterized. 
The  names  of  the  samples  were  given  according  to  the  gas  timing 
ratio  used  to  grow  them.  AES  measurements  show  that  the 
[0] :  [N  ]  ratio  of  the  samples  varied  with  the  gas  timing  used  (mono- 
tonically  but  not  linearly).  The  measured  XANES  spectra  and  the 
FEFF  simulations  are  shown  in  Fig.  7.  The  calculated  XANES  spectra 
of  wurtzite-InN  and  bixbyite  ln203  are  in  full  agreement  with  the 
measured  spectra  from  N30/00  and  NO/030  samples,  respectively. 
This  suggests  that  structurally  In  atoms  in  N30/00  are  mostly  four¬ 
fold  (as  those  in  InN)  and  in  NO/030  are  mostly  sixfold  (as  those  in 
ln203).  If  we  observe  the  change  in  the  spectra  while  the  oxygen 
content  is  increased  (from  the  bottom  curve  to  top  curve),  we  can 
see  that  all  five  features  marked  by  the  arrows  of  the  measured 
spectra  are  progressively  evolved.  However,  the  N30/05  spectrum, 
which  has  43%  O  content  [  1 0],  is  still  almost  identical  to  that  of  N30/ 
00  except  the  first  shoulder,  which  is  slightly  reduced.  This  indi¬ 
cates  that  most  of  In  atoms  in  this  sample  remained  fourfold  with 
some  0  atoms  substituting  on  the  N  sites.  The  simulated  spectrum 
of  InN0.600.4  alloy  in  wurtzite  structure  was  calculated  using  FEFF 
codes  to  test  our  assumption.  The  calculated  spectrum  is  almost 
overlapped  with  that  of  (calculated)  pure  InN  except  the  first  shoul¬ 
der  is  reduced;  in  full  agreement  with  experimental  observation 
and  DFT  PDOS  discussed  in  previous  section.  As  explained  in  the 
previous  section,  the  reduction  of  the  first  shoulder  occurs  because 
parts  of  the  CB  are  occupied.  For  samples  with  higher  0  contents 
(N30/010  and  N30/020),  the  spectra  show  a  mixed  signature  be¬ 


Photon  Energy  (eV) 


tween  those  of  wurtzite-InN  and  bixbyite  ln203.  This  means  that 
there  are  both  fourfold  and  sixfold  In  atoms  in  these  samples.  The 
excess  0  atoms  above  the  substitutional  solubility  limit  (assumed 
to  be  ~40%)  of  InN  can  form  a  phase  separated  ln203-like  structure 
that  lead  to  a  formation  of  sixfold  In  atoms.  Based  on  this  model,  the 
composition  in  the  sample  can  be  written  as  (InNo.eOo^MInCh .5)i_f, 
where  the  F  values  are  0.49  and  0.38  for  N30/010  and  N30/020, 
respectively  [10].  The  simulated  spectra  calculated  using  weight- 
averaged  between  InN06Oo.4  and  ln203  are  shown  in  Fig.  7.  The  sim¬ 
ulated  spectra  based  on  this  model  are  consistent  with  measured 
spectra  (N30/010  and  N30/020). 

For  ln203,  since  there  are  two  species  of  nonequivalent  In 
atoms,  the  simulated  spectrum  shown  in  Fig.  7b  (top  curve)  is  ob¬ 
tained  from  the  weight-averaged  (Inavg)  between  the  spectra  of 
both  In  species  (Ini  and  In2).  As  can  be  seen  in  Fig.  la,  the  local 
structures  are  quite  different  (although  both  Ini  and  In2  are  six- 
folded),  one  might  be  curious  how  similar  the  spectra  of  Ini  and 
In2  are.  The  spectra  obtained  from  the  two  species  of  In  atoms  in 
ln203  are  quite  similar,  as  shown  in  Fig.  8b.  This  further  supports 
the  fact  that  the  shape  of  the  spectra  depends  strongly  on  the  coor¬ 
dination  number  (fourfold  or  sixfold)  but  only  weakly  on  the  de¬ 
tailed  geometry  distortion.  For  further  test,  the  effects  of  the 
substitution  of  a  neighboring  N  by  0  (or  0  by  N)  were  studied.  First, 
we  replaced  one  N  atom  surrounding  an  In  atom  with  an  0  atom  in 
a  wurtzite-InN  supercell  (72-atom).  Second,  all  four  N  atoms  sur¬ 
rounding  the  In  atom  were  replaced  by  0  atoms  in  a  wurtzite- 
InN  supercell.  The  simulated  In  L3-edge  XANES  spectra  of  the  In 
atoms  that  their  neighbors  were  replaced  are  shown  (labeled  10N 
and  404N  for  the  cases  that  one  and  four  neighbors  are  replaced, 
respectively)  in  comparison  with  that  of  pure  InN  in  Fig.  8a.  It  is 
clear  that  even  all  four  N  neighbors  are  replaced,  the  features  of 
the  spectrum  still  remain  the  same  as  those  of  pure  InN.  We  have 
also  tested  the  bixbyite  (sixfold  In)  structure.  We  replaced  one  oxy¬ 
gen  atom  by  a  nitrogen  atom  in  the  ln203  bixbyite  structures.  The 
simulated  XANES  spectrum  of  the  In  atom  that  its  neighboring  0 
was  replaced  by  N,  is  shown  (labeled  1N0)  in  comparison  with  that 
of  pure  ln203  in  Fig.  8b.  It  is  clear  that  the  XANES  features  are  not 
very  sensitive  to  anion  substitutions  but  very  sensitive  to  the  coor¬ 
dination  number  of  In  (fourfold  or  sixfold). 


Fig.  7.  The  measured  and  simulated  In  L3-e dge  XANES  spectra  of  all  five  samples. 
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Fig.  8.  (a)  The  simulated  In  L3-edge  XANES  spectra  of  In  atoms  in  pure  wurtzite  InN  (bottom),  in  wurtzite  InN  with  one  N  replaced  by  0  (middled),  and  in  wurtzite  InN  with 
four  N  replaced  by  0  (top),  (b)  The  spectra  In  atoms  in  pure  bixbyite  ln203  (bottom  three  curves)  and  in  bixbyite  ln203  with  one  0  replaced  by  N  (top). 


4.  Conclusion 

Synchrotron  X-ray  absorption  near  edge  structures  of  In  L3-edge 
were  used  to  study  the  local  geometry  of  indium  atoms  in  indium 
oxynitride  thin  films  with  varied  O  contents.  Calculations  showed 
that  the  spectra  features  are  sensitive  to  the  coordination  number 
of  the  In  atom  (fourfold  for  wurtzite  structure  or  sixfold  for  bixby¬ 
ite  structure)  but  not  very  sensitive  to  the  species  (0  or  N)  of  the 
first  neighbors  such  that  interchanging  N  with  O  leads  to  almost 
unnoticeable  changes  in  the  XANES  spectra.  The  solubility  limit 
of  oxygen  in  wurtzite  InN  is  estimated  to  be  ~40%  before  the  crys¬ 
tal  structure  drastically  changed.  The  band  structures  and  density 
of  states  of  InN,  ln203,  and  their  alloys  are  shown  and  discussed. 
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Chen  et  al.  [Phys.  Rev.  Lett.  96  (2006)  035508]  experimentally  observed  vibrational  signatures  related  to 
defects  in  oxygen-doped  CdTe  using  ultrahigh  resolution  Fourier  transform  infrared  (FTIR)  spectroscopy. 
They  observed  an  absorption  peak  at  350  cm  ’.In  addition,  for  samples  grown  under  certain  conditions, 
they  observed  two  higher  frequency  peaks  (1097  and  1108  cm  ’)  at  low  temperature  that  merged  into 
one  at  room-temperature.  They  attributed  the  low-frequency  peak  (350  cm  ’)  to  the  vibration  of  0Te 
and  the  two  higher  frequency  peaks  to  the  vibrational  modes  of  a  Ofe-Vcd  complex.  Subsequently,  they 
reported  similar  modes  around  1100  cm  1  in  O-doped  CdSe  [Phys.  Rev.  Lett.  101  (2008)  195502]  which 
were  attributed  to  an  0Se-VCd  complex.  We  employed  first-principles  DFT  calculations  to  calculate  the 
vibrational  modes  of  0Te  and  0Te-VCd  complex  in  CdTe.  Our  calculations  show  that  the  350  cm  1  mode 
is  consistent  with  0Te.  However,  the  frequencies  of  the  modes  around  1100  cm  1  are  more  than  twice 
the  expected  frequencies  for  0Te-Vcd  complexes  in  CdTe  (or  0Se-Vcd  in  CdSe),  indicating  that  the  0Te- 
VCd  complex  cannot  be  the  cause  of  the  observed  1100  cm  1  modes.  A  search  for  a  new  defect  model 
is  in  order. 

©  2010  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

CdTe  is  recognized  as  a  promising  semiconductor  because  it  can 
be  doped  both  n-type  and  p-type.  The  role  of  defects  in  CdTe  has 
been  theoretically  studied  by  Wei  and  Zhang  [1].  CdTe  and  its  al¬ 
loys  have  been  used  for  applications  such  as  room-temperature 
X-ray  and  y-ray  detectors  [2,3],  The  Bridgman  technique  [4,5] 
can  be  used  to  grow  high  quality  CdTe  bulk  crystal  samples. 

Recently,  local  vibrational  modes  (LVM)  related  to  oxygen 
impurities  in  CdTe  have  been  experimentally  studied  using  Fourier 
transform  infrared  (FTIR)  spectroscopy  [4,6].  Depending  on  the 
growth  and  doping  conditions,  low  temperature  FTIR  measure¬ 
ments  reveal  three  LVMs  related  to  oxygen:  one  low-frequency 
mode  at  350cnrr’  and  two  high-frequency  modes  at  v !  =  1097 
and  v2  =  1108  cur1.  As  the  measuring  temperature  is  increased, 
the  two  high-frequency  modes  (v!  and  v2)  merge  into  one  at 
room-temperature.  Chen  et  al.  [4,6]  assigned  the  low-frequency 
mode  to  an  oxygen  substituting  on  the  Te  site  (0Te),  and  the  two 
high-frequency  modes  to  a  complex  formed  by  an  oxygen  substi¬ 
tuting  on  a  Te  site  and  a  neighboring  Cd-vacancy  (0Te-VCd)  [4].  La¬ 
ter,  they  also  found  similar  high-frequency  modes  around 
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1100  cur1  in  O-doped  CdSe  [7]  that  they  attributed  to  an  0Se- 
Va  complex.  Chen  et  al.  believe  that  when  forming  a  complex  with 
Vca,  the  frequency  of  0Te  can  increase  by  almost  a  factor  of  three. 
They  attributed  the  increase  in  the  frequency  to  an  increase  in 
the  bond  strength  between  the  O  atom  and  its  neighbors  as  the 
number  of  Cd  neighbors  is  reduced  from  four  (for  0Te)  to  three 
(for  0Te-Vcd). 

To  understand  these  vibrational  modes,  we  have  investigated 
O-related  centers  in  CdTe  using  density  functional  calculations. 
The  LVM  frequencies  of  0Te  and  0Te-VCd  defects  were  calculated 
by  using  force-response  matrix  (dynamical  matrix)  calculations. 
The  harmonic  approximation  is  used.  We  will  show  that  the  calcu¬ 
lated  LVM  of  0Te  is  consistent  with  the  experimental  value.  How¬ 
ever,  the  calculated  LVM  frequencies  of  0Te-VCd  are  much  lower 
than  the  Vx  and  v2  modes  observed  by  Chen  et  al.  [4],  Although 
the  O-Cd  bonds  are  slightly  stronger,  the  frequency  increases  by 
only  38%,  not  by  a  factor  of  three.  Therefore,  it  is  clear  that  the 
0Te-VCd  complex  is  not  the  correct  model  to  explain  the  observed 
LVMs.  Preliminary  results  of  our  work  have  been  published  as  a 
comment  [5[.  Here,  we  present  the  computational  method  and  de¬ 
tails  of  the  LVMs.  The  two  high  frequencies  observed  by  Chen  et  al. 
[4]  should  be  attributed  to  other  defects  related  to  the  introduction 
of  O  in  CdTe,  not  to  0Te-VCd-  We  note  that  the  high  value  of  the  fre¬ 
quency  makes  hydrogen  (which  could  be  related  to  the  introduc¬ 
tion  of  O  during  growth)  a  prime  candidate. 
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2.  Calculations  of  local  vibrational  mode 

Our  work  is  based  on  density  functional  theory  within  the  local 
density  approximations  (LDA).  For  the  electron-ion  interactions, 
we  used  projector  augmented  wave  potentials  [8],  as  implemented 
in  the  VASP  code  [9,10],  with  an  energy  cutoff  of  500  eV  for  the 
plane-wave  basis  set.  The  calculated  lattice  constant  of  bulk  CdTe 
is  6.421  A,  in  good  agreement  with  the  experimental  value  of 
6.477  A  [1 1  ].  To  calculate  the  equilibrium  lattice  constant,  calcula¬ 
tions  of  bulk  CdTe  at  various  lattice  constants  were  performed,  and 
fitting  to  an  equation  of  state  [12]  identifies  the  lowest  energy 
structure. 

For  defects,  a  supercell  approach  is  used  [13].  In  principle,  to 
simulate  an  isolated  defect,  the  larger  supercell  size,  the  better. 
Previously  we  have  shown  that  a  32-atom  is  sufficient  to  suppress 
the  interaction  between  the  defect  and  its  images  [14].  Here,  we 
used  a  64-atom  supercell,  which  is  a  2  x  2  x  2  repetition  of  the 
CdTe  conventional  8-atom  unit  cell.  All  atoms  in  the  supercell  were 
allowed  to  relax  until  the  Hellmann-Feynman  forces  [15]  become 
less  than  1CT3  eV/A.  The  Monkhorst-Pack  scheme  [16]  with  a  sam¬ 
pling  mesh  of  2x2x2  special  k-points  is  used  for  fc-space 
integrations. 

The  LVMs  of  0  defects  were  calculated  using  the  so-called  fro¬ 
zen  phonon  approach  [17],  First,  for  each  defect  configuration, 
the  structure  is  allowed  to  fully  relax.  Then  we  slightly  shifted  each 
and  every  atom  in  the  supercell  along  all  three  axes,  one  at  a  time, 
and  the  dynamical  matrix  can  be  constructed  from  the  forces  on  all 
atoms  in  response  to  these  displacements.  Within  the  harmonic 
approximation,  the  total  energy  of  a  supercell  with  small  displace¬ 
ments  of  atoms  from  their  equilibrium  positions  can  be  written 
as 

«  =  a) 

ij>,0 

where  dx(i )  is  the  displacement  of  atom  i  from  its  equilibrium  posi¬ 
tion  in  the  direction  a  and  <Pa/5(i,  j)  is  the  real-space  force  constant 
matrix.  A  small  displacement  d/;(j)  of  atom  j  in  direction  fl  induces 
a  force  on  atom  i  in  direction  a  described  by 

j)df(j).  (2) 

Based  on  Eq.  (2),  we  can  construct  the  complete  real-space  force 
constant  matrix  c Pa/i{uj).  This  is  done  by  creating  a  small  displace¬ 
ment1 *  d0  (we  use  dp(j)  =  d0  =  0.01  x  a  where  a  is  the  lattice  constant) 
of  each  atom  j  in  the  supercell  in  three  orthogonal  directions, 
p  =  1,  2, 3,  one  atom  and  one  direction  at  a  time.  For  each  calculation, 
the  (Hellmann-Feynman)  forces  on  all  atoms  in  all  three  directions 
Fa(i)  are  recorded  and  the  components  of  the  force  constant  matrix 
can  be  obtained  by  using  Eq.  (2).  To  reduce  from  the  impact  of  anhar- 
monic  effects,  each  component  of  the  force  constant  matrix  is  ob¬ 
tained  from  an  average  of  two  calculations:  one  with  +d0 
displacement  and  another  with  -d0  displacement.  The  dynamical 
matrix  can  be  easily  calculated  from  the  force  matrix  according  to 
the  equation, 

=  (MdW,r1/2<MU),  (3) 

where  is  the  mass  of  atom  i.  The  diagonalization  of  the  dynam¬ 
ical  matrix  yields  the  eigenvalues  ek  =  col,  and  eigenvectors  \uk). 
The  eigenvalues  are  the  frequencies  of  the  vibrational  modes.  The 


1  The  displacement  should  be  large  enough  to  allow  an  accurate  determination  of 
the  energy  increase  (or  response  forces)  in  the  calculations.  However,  it  should  be 

small  enough  to  stay  within  the  range  of  realistic  vibrational  amplitudes  where  the 

harmonic  approximation  is  expected  to  be  valid. 


eigenvector  of  each  mode  describes  the  displacements  of  each  atom 
in  the  system. 


3.  Results  and  discussion 

Fig.  la  shows  the  crystal  structure  of  CdTe.  The  0Te  center  is  cal¬ 
culated  by  replacing  one  Te  atom  in  the  64-atom  CdTe  supercell 
with  an  O  atom  (0Te)  and  allowing  full  relaxation.  We  found  that 
the  O  atom  forms  bonds  with  its  four  Cd  neighbors  with  an  equal 
distance  (Fig.  lb).  The  bond  distance  of  O-Cd  is  shorter  than  that 
of  Cd-Te  bond  in  bulk  CdTe  by  19%. 

The  complex  between  0Te  and  a  Cd-vacancy  (Che-Vcd)  is  created 
by  removing  a  Cd  neighbor  of  0Te.  Under  n-type  and  semi-insulat- 
ing  conditions,  we  found  that  the  0Te-Vcd  complex  is  stable  in  the 
2~  charge  state.  Because  CdTe  samples  used  in  experiments  are 
generally  semi-insulating  [4],  we  reported  only  the  2~  charge  state 
in  this  paper  (Note  that  other  two  possible  charge  states,  i.e.  neu¬ 
tral  and  1~,  have  very  similar  bond  structures  and  vibration  fre¬ 
quencies.).  Since  one  Cd-0  bond  is  missing,  the  O  atom  moves 
closer  to  the  three  remaining  Cd  neighbors.  As  a  result  the  O  atom 
is  moved  in  the  direction  away  from  the  center  of  the  vacancy  to¬ 
ward  the  center  of  the  three  remaining  Cd  neighbors,  leading  to 
shorter  bond  distances  between  the  O  atom  and  the  three  remain¬ 
ing  Cd  neighbors.  The  remaining  three  Cd-0  bonds  are  equal  in  dis¬ 
tance  and  32%  shorter  than  the  Cd-Te  distance  in  bulk  CdTe 
(Fig.  lc).  Note  that  the  O  atom  is  equidistant  from  the  three  Cd 
neighbors  and  only  0.26  A  away  from  the  plane  formed  by  those 
Cd  atoms. 

The  LVM  frequencies  of  0Tc  and  the  0Te-Vcd  complex  were  cal¬ 
culated  by  using  the  full  dynamical  matrix  calculations  for  a  64- 
atom  supercell  as  described  in  the  previous  section.  To  identify 
the  modes  localized  on  the  O  impurity,  we  calculated  the  sum- 
square  of  the  three  eigenvector  components  associated  with  the 
O  atom  (Jjatom=ou2  or.  in  short,  J2uo)-  Note  that  the  sum-square 
of  all  components  is  normalized  to  1.  The  for  all  modes  are 
plotted  as  histograms  in  Fig.  2.  For  0Te,  we  can  clearly  see  that 
the  triply  degenerate  mode  at  338  cm-1  is  localized  on  the  O  atom. 
The  mode  is  triply  degenerate  because  the  O  atom  has  fourfold 
symmetry.  Under  a  simple  harmonic  approximation,  the  vibrations 
of  the  O  atom  in  all  three  axes  are  therefore  equivalent.  An  illustra¬ 
tion  of  the  set  of  eigenvectors  describing  these  degenerate  modes 
is  shown  in  Fig.  3a.  The  value  of  J2  Ug  for  this  mode  is  0.95,  indicat¬ 
ing  a  very  strong  localization  of  the  mode.  The  calculated  LVM  fre¬ 
quency  of  338  cnr1  is  in  good  agreement  with  the  observed  value 
of  350  cnr1  [4[.  For  the  0Te-VCd  complex  (Fig.  2b),  we  find  a  doubly 
degenerate  mode  (at  467  cnr1)  and  a  singlet  mode  (at  197  cnr1) 
that  are  localized  on  the  O  atom.  The  values  of  J2  Uo  are  0-94  and 
0.34  for  the  doublet  and  singlet  states,  respectively.  These  modes 
are  derived  from  the  triply  degenerate  mode  of  0Te.  The  introduc¬ 
tion  of  VCd  leads  to  a  symmetry  breaking  that  separates  the  modes 
into  a  doublet  and  a  singlet.  The  doublet  (467  cm-1)  is  associated 
with  the  in-plane  vibration;  the  O  motion  is  illustrated  in  Fig.  3b. 
The  singlet  (197  cm-1)  is  associated  with  the  out-of-plane  vibra¬ 
tion  of  the  O  atom  accompanied  by  the  outward  breathing  of  the 
surrounding  Cd  neighbors  as  illustrated  in  Fig.  3c.  The  increase  in 
the  frequency  of  the  doublet  (in  comparison  to  that  of  0Te)  by 
38%  can  be  attributed  to  the  stronger  Cd-0  bond  in  the  complex. 
Similarly,  the  decrease  in  the  frequency  of  the  singlet  (in  compar¬ 
ison  to  that  of  0Te)  by  42%  can  be  attributed  to  the  missing  Cd-0 
bond  (in  the  direction  of  vibration)  of  the  complex. 

The  calculated  frequencies  of  197  cm-1  for  the  out-of-plane 
vibration  and  467  cnr1  for  the  in-plane  vibration  are  completely 
inconsistent  with  the  two  high  LVM  frequencies  of  1097  and 
1108  cnr1  assigned  by  Chen  et  al.  [4]  for  out-of-plane  and  in-plane 
vibration  of  the  0Te-Vcd  complex.  While  the  presence  of  a  Cd-va- 
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Fig.  1.  The  local  structures  of  (a)  bulk  CdTe,  (b)  0Te,  and  (c)  0Te-Vcd  complex  in  CdTe.  The  bond  distances  are  given  as  percentage  differences  from  the  Cd-Te  bond  distance 
(calculated  value  =  2.781  A)  in  bulk  CdTe. 
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Fig.  2.  The  sum-square  of  the  three  eigenvector  components  associated  with  the  0  atom  (Satom.ou2)  for  vibrational  modes  for  (a)  0Te  and  (b)  the  Olr.-V,(i  complex  in  a  64- 
atom  CdTe  supercell. 


Fig.  3.  Schematic  illustrations  of  the  vibration  vectors  describing  the  local  vibrational  modes  associated  with  O  in  0Te  and  in  the  Ole-14ci  complex  in  CdTe.  (a)  The  triply 
degenerate  mode  of  0Te.  (b)  The  doubly  degenerate  mode  and  (c)  the  singlet  mode  of  the  0Te-Vcd  complex  in  CdTe.  In  (a)  and  (b)  each  arrow  represents  the  vibration  vector  of 
the  O  atom  in  each  of  the  degenerate  mode.  Components  on  other  atoms  are  small  and  not  shown.  In  (c),  the  four  arrows  show  the  vibration  vectors  of  the  O  atom  and  its  Cd 
neighbors. 


cancy  near  the  substitutional  oxygen  indeed  leads  to  an  increase  in 
the  vibrational  frequency  of  the  in-plane  vibrational  mode,  the 
magnitude  of  this  increase  is  far  smaller  than  Chen  et  al.  assumed. 
The  frequency  is  increased  from  338  cnr1  to  467  cnr1,  a  mere 
38%,  very  different  from  the  228%  assumed  in  Chen  et  al.’s  attribu¬ 
tion.  Moreover,  the  out-of-plane  frequency  does  not  increase  at  all; 
instead,  it  is  lowered  to  197  cm-1.  Based  on  our  calculations,  we 
can  therefore  safely  conclude  that  the  high-frequency  modes 


(1097  and  1108  cnrr1)  observed  in  the  experiment  of  Chen  et  al. 
[4]  are  not  associated  with  a  0Te-Vcd  complex. 

More  recently,  similar  LVM  frequencies  (1094,  1107,  and 
1126  cnr1)  were  observed  in  oxygen-doped  CdSe  by  the  same 
experimental  group  [7].  They  similarly  assigned  these  LVM  fre¬ 
quencies  to  an  Ose-Vcd  complex  and  explained  that  an  additional 
mode  (they  observed  three  high-frequency  modes  instead  of  the 
two  modes  in  the  case  of  CdTe)  is  the  result  of  the  lower  crystal 
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symmetry  of  CdSe.  For  completeness,  we  also  calculated  the  vibra¬ 
tional  frequencies  of  0Se  and  the  0Se-Vcd  complex  in  CdSe  using  the 
same  approach  as  for  the  O-related  centers  in  CdTe.  We  obtained 
the  LVM  frequency  of  0Se  at  370  cm1,  comparable  to  (and  slightly 
higher  than)  that  of  0Te  in  CdTe.  For  the  0Se-Vcd  complex,  we  ob¬ 
tained  LVM  frequencies  of  512,  507,  and  223  cirr1  again  compara¬ 
ble  to  (and  consistently  higher  than)  the  corresponding  calculated 
frequencies  for  the  0|e-Vtd  complex  in  CdTe.  Because  the  calcu¬ 
lated  frequencies  are  inconsistent  with  the  observed  values  of 
1094,  1107,  and  1126  cm~\  we  can  again  safely  conclude  that 
the  observed  modes  are  not  associated  with  Ose-Vcd  complex  in 
CdSe. 

We  acknowledge  that  first-principles  calculations  of  vibrational 
frequencies  have  error  bars,  but  in  general  the  calculated  frequen¬ 
cies  agree  very  well  with  experimental  values  [18-23],  Deviations 
from  experimental  values  by  more  than  15%  (or  about  150  cm-1) 
are  very  rare.  The  huge  qualitative  difference  between  the  frequen¬ 
cies  for  the  Cd-vacancy  related  complexes  and  the  modes  observed 
by  Chen  et  al.  [4,7]  therefore  clearly  indicates  that  the  attribution 
of  these  modes  to  0Te(orSe)-^cd  complexes  cannot  be  correct. 


4.  Conclusion 

Based  on  first-principles  density  functional  calculations,  the 
vibrational  frequencies  associated  with  0Te  and  0Te-Vcd  complexes 
in  CdTe  and  0Se  and  0Se-Vcd  complexes  CdSe  were  calculated  and 
used  to  interpret  recent  experimental  measurements  by  Chen  et  al. 
[4,6].  In  the  experiment,  a  low-frequency  local  vibrational  mode 
(350  cirr1)  and  two  high-frequency  modes  (1097  and 
1108  cirr1)  are  observed  in  oxygen-doped  CdTe  samples.  The 
low-frequency  mode  has  been  assigned  to  0Te  and  the  two  high- 
frequency  modes  have  been  assigned  to  a  0Te-Vcd  complex.  In  a 
more  recent  experiment  on  CdSe  by  the  same  group,  similar  modes 
are  observed  and  assigned  to  0Se  and  Ose-Vcd.  Our  calculated  fre¬ 
quency  of  0Te  in  CdTe  is  338  cur1  which  is  in  a  reasonable  agree¬ 
ment  with  the  measured  value  of  350citt'.  However,  the 
calculated  frequencies  of  the  Oxe-Vca  complex  are  197  and 
467  cm~\  far  lower  than  the  observed  values  (1097  and 
1108  cirr1)  and  strongly  indicating  that  the  assignment  of  the 
high-frequency  modes  to  the  0Te-Vcd  complex  is  incorrect.  The 
1097  and  1108  cm~'  modes  must  come  from  other  defects,  possi¬ 


bly  hydrogen  related.  Similar  conclusions  apply  to  the  case  of  oxy¬ 
gen  in  CdSe. 
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Direct  enumeration  method  is  an  efficient  way  for  scanning  alloy  properties  by  computing  all  possible 
configurations.  This  method  thoroughly  covers  all  possible  configurations  without  bias  and  provides  an 
informative  trend  of  alloy  properties  on  ordering  patterns.  In  an  actual  study,  the  number  of  possible  con¬ 
figurations  increases  rapidly  with  the  size  of  the  unit  cell  and  usually  a  reasonable  size  that  give  a  suffi¬ 
ciently  large  number  of  configurations  are  chosen.  In  this  work,  the  convergence  of  the  bandgap  range 
with  respect  to  the  unit  cell  size  of  AlGalnP  alloy  is  studied  up  to  8  cation  atoms  per  unit  cell.  It  is  found 
that  the  bandgap  range  already  converges  to  within  0.1  eV  when  the  unit  cell  size  is  4  cation  atoms.  Inter¬ 
estingly,  we  also  found  that  the  lowest  bandgap  value  of  GalnP  alloy  is  achieved  already  in  a  small  cell  (2 
cations  in  the  unit  cell).  The  cause  for  this  special  small  bandgap  case  is  discussed. 

©  2010  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

Alloying  provides  a  flexible  way  to  engineer  electronics  and 
optical  properties  of  semiconductors  to  suit  the  need  of  new  gen¬ 
eration  device  applications  [1-3].  The  properties  of  semiconductor 
alloys,  especially  those  of  III— V  alloys,  have  been  extensively  stud¬ 
ied  in  the  last  decade  [1-4].  Several  schemes  were  proposed  to  ex¬ 
plain  the  relationship  between  the  properties  and  the  alloy 
composition.  The  simplest  and  most  known  scheme  (the  Vegard’s 
rule)  assumes  a  linear  dependence  between  the  interested  prop¬ 
erty  and  the  alloy  composition  [5,6].  In  reality,  most  properties 
are  not  linearly  dependent  with  the  alloy  composition.  As  a  result, 
non-linear  terms  have  been  proposed  to  quantify  the  deviation 
from  the  Vegard’s  trend.  The  most  known  improved  scheme  added 
a  bowing  parameter  [7[.  However,  neither  the  Vegard's  rule  nor  its 
variants  takes  into  account  the  ordering  of  the  constituent  atoms  in 
the  alloy.  It  is  important  to  note  that,  for  a  given  alloy  composition, 
there  are  virtually  infinite  possible  arrangements  of  the  atoms  in 
the  alloy.  From  the  theoretical  point  of  view,  the  direct  enumera¬ 
tion  method  (DEM)  has  been  used  to  study  alloy  systems  [8-11], 
The  important  elements  of  the  method  are  the  algorithm  to  gener¬ 
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ate  alloy  configurations  [10,12,13]  and  the  practical  method  to  cal¬ 
culate  their  properties  [8,10,11,14-17].  The  algorithm  to  generate 
alloy  configurations  for  binary  alloy  with  face-centered  cubic 
(FCC)  and  body-centered  cubic  (BCC)  lattice  have  been  proposed 
by  Ferreira  et  al.  [10].  Later,  an  efficient  algorithm  for  generating 
alloy  configurations  has  been  introduced  by  Hart  and  Forcade 
[18]  where  the  run  time  scales  linearly  with  the  number  of  distinct 
configurations.  Recently,  an  algorithm  for  generating  alloy  config¬ 
urations  of  a  non-primitive  lattice  and  application  to  hexagonal 
close  packed  alloys  have  been  developed  [19],  The  DEM  has  been 
used  to  study  the  bandgap  properties  calculated  by  empirical 
pseudopotential  method  (EPM)  [14-17,20-23]  of  several  classes 
of  semiconductor  alloys;  including  the  ternary  [9-11]  and  quater¬ 
nary  [24]  alloys.  An  advantage  of  DEM  is  that  it  relates  the  inter¬ 
ested  property  to  the  alloy  configuration  (atomic  ordering)  in 
additional  to  the  composition.  From  EPM  calculation,  we  have  pre¬ 
viously  shown  that  the  bandgap  of  AlGalnP  alloy  depends  strongly 
on  the  atomic  ordering  in  addition  to  the  alloy  composition  [24].  In 
that  work,  we  have  also  shown  that  most  of  the  alloy  configura¬ 
tions  have  the  bandgaps  that  are  smaller  than  the  values  predicted 
by  the  Vegard’s  rule.  The  ordering  causes  a  reduction  in  the  band- 
gap;  pushing  it  below  the  Vegard’s  trend  [24-26].  The  lowering  of 
the  bandgap  due  to  atomic  arrangements  emphasizes  the  role  of 
the  alloy  configurations  (in  addition  to  the  alloy  composition)  in 
bandgap  engineering.  In  this  work  and  our  previous  work  [24], 
we  focused  our  attentions  in  studying  a  complete  set  of  theoretical 
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possible  arrangements  within  a  given  cell  size.  These  complete  re¬ 
sults  allow  us  to  investigate  how  the  band  properties  relate  to  the 
atomic  arrangements.  However,  we  should  keep  in  mind  that  many 
of  the  configurations  studied  might  have  a  very  high  energy  and 
difficult  to  fabricate. 

Because  there  are  many  distinct  alloy  configurations  for  a  given 
alloy  composition,  DEM  gives  the  possible  range  of  properties  (for 
e.g.,  the  bandgap  and  the  electron  effective  mass)  at  a  given  alloy 
composition,  as  well  as  the  relationship  between  the  ordering 
and  a  given  property.  While  a  larger  cell  size  gives  more  possible 
configurations  (allowing  a  more  complete  study)  and  finer  compo¬ 
sition  steps,  the  computational  demand  is  increased  rapidly.  There¬ 
fore,  it  is  important  to  investigate  the  size  of  unit  cell  that 
sufficiently  gives  the  converged  results.  In  this  work,  we  focus  on 
the  quaternary  alloys  AlxGayIni_x_yP  (also  called  pseudoternary  al¬ 
loys  that  contain  three  mixed-cations  in  the  FCC  sublattice). 

2.  Calculation 

The  alloy  configurations  of  AlGalnP  are  generated  by  applying 
DEM  [9,11,24]  to  a  zincblende-based  alloy  system.  The  possible 
number  of  alloy  configurations  for  a  unit  cell  containing  N  cation 

atoms  is  given  by  the  multichoose  formula  [27],  +  ^  'j  — 

(n  -  1,N)!,  with  n  being  the  number  of  elements.  For  AlGalnP  alloy, 
n  =  3.  A  subroutine  in  the  Alloy  Theoretic  Automated  Toolkit 
(ATAT)  package  [28,29]  was  employed  to  systematically  generate 
the  configurations.  The  detail  of  alloy  configurations  generated 
by  ATAT  is  described  elsewhere  [24].  Note  that,  some  of  the  gener¬ 
ated  configurations  are  equivalent  by  symmetries  and  the  redun¬ 
dant  configurations  were  removed  to  reduce  computation 
demand.  Using  the  ATAT  codes,  the  cell  size  of  up  to  8  cation  atoms 
per  unit  cell  (Nmax  =  8)  provides  a  total  of  9808  distinct  configura¬ 
tions.  They  are  composed  of  1887  pseudobinary  alloy  configura¬ 
tions  from  the  three  parent  binary  compounds  (3  x  629  AlGaP, 
GalnP  and  AllnP)  and  7918  pseudoternary  alloy  configurations 
[see  Ref.  [24]  (Table  1)  for  detail]. 

Because  the  composition  of  the  AlxGayIni_x_yP  alloys  is  defined 
by  two  variables  (x  and  y),  the  composition  space  can  be  presented 
in  a  two-dimensional  space.  We  have  previously  shown  that  the 
two-dimension  composition  space  of  this  ternary  alloy  system  is 
well  described  by  using  the  transformed  coordinates,  u=x+y/2 
and  v  =  V3y/2  [24].  This  leads  to  a  regular  triangle,  as  shown  in 
Fig.  1 .  The  parent  compounds  are  located  at  the  vertices  and  the 
circles  represent  the  alloy  compositions  that  can  be  obtained  by 
each  unit  cell  size. 

To  reduce  the  computation  demand,  alloy  properties  of  all  gen¬ 
erated  configurations  are  calculated  by  using  EPM  [14,15].  The 
structural  relaxation  was  done  by  the  valence  force  field  method 
(VFF)  [30-33],  In  EPM,  the  electronic  structures  are  obtained  by 
solving  the  simplified  Schrodinger  equation 


Fig.  1.  The  composition  space  of  the  pseudoternary  alloy  system.  The  colored 
circles  show  the  possible  compositions  given  by  each  unit  cell  size  (i.e. 
N  =  2, 3,4,  5, 6, 7,  and  8).  Circles  on  the  three  edges  of  the  triangle  correspond  to 
the  pseudobinary  alloys. 

V2  +  £  -  Ra,„)W)  =  £iUr)  (1) 

Z  «,n  J 

where  the  effective  pseudopotential  is  summed  over  all  individual 
atoms  n  (a  denotes  atomic  type).  Va  is  the  screened  potential  for 
atomic  type  a  and  is  the  atomic  position  of  atom  n  (type  a)  ob¬ 
tained  after  VFF  relaxations.  One  advantage  of  EPM  (over  first-prin¬ 
ciples  density  functional  calculations)  is  that  it  provides  the 
bandgaps  that  do  not  need  corrections.  This  is  because  the  parame¬ 
ters  are  already  fitted  to  the  experimental  bandgaps.  The  cutoff  en¬ 
ergy  of  the  generated  pseudopotentials  is  68  eV.  The  details  of  the 
fitting  parameters  for  the  atomic  potentials  of  A1P,  GaP,  and  InP  com¬ 
pounds  can  be  found  in  Ref.  [16].  Details  of  EPM  method  [14-16]  and 
electronic  structure  computation  are  described  elsewhere  [24]. 

For  each  set  of  the  configurations  in  a  unit  cell  size,  the  maxi¬ 
mum  and  the  minimum  bandgaps  are  determined.  The  possible 
range  of  bandgaps  is  the  difference  between  the  maximum  and 
the  minimum  values  at  a  given  composition.  Because  the  ordering 
of  cation  atoms  can  cause  a  bandgap  reduction  [24],  the  larger  cell 
size,  which  facilitates  more  degree  of  ordering,  is  likely  to  provide  a 
smaller  minimum  bandgap  than  the  smaller  cell  size.  However,  it  is 
expected  that  the  minimum  value  of  bandgap  would  converge  at  a 
sufficiently  large  cell.  In  this  work,  the  largest  cell  size  used  is  the 
unit  cell  with  8  cations.  The  bandgap  minima  are  calculated  for  cell 
size  of,  JVmax  =  4,  5,  6,  7,  and  8.  The  upper  bound  of  bandgap  is  set  by 
the  Vegard’s  rule  which  represent  the  disordered  alloys  [24]. 
According  to  Vegard’s  rule,  the  (upper  bound)  bandgap  of  Alx_ 
GayIni_x_yP  is 

£gap(x,y)  =  x£gap(AlP)  +y£gap(GaP)  +  (1  -  x  -y)Egap(lnP)  (2) 


Table  1 

The  selected  alloy  configurations  with  ordering  in  the  [111]  directions.  Only  the  positions  of  the  cations  are  listed  (the  position  of  P  atoms  are  shifted  by  (0.25, 0.25, 0.25)  from 
the  cation  positions).  The  lattice  vectors  and  the  cation  positions  are  in  the  Cartesian  coordinates  in  the  unit  of  lattice  constant. 


Alloy  configuration 

Lattice  vector 

Cation  position 

Alloy  ordering 

I 

(-0.5,  0.5, 1)  (-0.5, 1, 0.5)  (-1,  0.5,  0.5) 

Ga:  (-2,2,2) 

In:  (-1,1,1) 

(GaPhIKInP), 

II 

(-1.5, 1, 1.5)  (-1, 1.5, 1.5)  (-1.5, 1.5, 1) 

Ga:  (-1,1,1),  (-4, 4,4) 

In:  (-2,2,2),  (-3,  3,3) 

(GaP)2||(InP)2 

Ilia 

(-0.5,  0,  -0.5)  (0,  -0.5,  0.5)  (-2,  2,  2) 

Ga:  (-0.5,  0,  0.5),  (-1,0.5,  0.5),  (-2.5, 1.5,2) 

In:  (-1.5,  0.5, 1),  (-1.5, 1, 1.5),  (-2, 1.5, 1.5) 

(GaP)3||(InP)3 

mb 

(-0.5,  0,  -0.5)  (0,  -0.5,  0.5)  (-2,  2,  2) 

Ga:  (-0.5,  0,  0.5),  (-1.5,  0.5, 1),  (-2.5, 1.5,  2) 

In:  (-1, 0.5,  0.5),  (-1.5, 1, 1.5),  (-2, 1.5, 1.5) 

(GaP)1||(InP)1||(GaP)2||(InP)2 

IV 

(-2.5,  2.5,  3)  (-2.5, 3, 2.5)  (-3, 2.5,  2.5) 

Ga:  (-1,1,1),  (-3, 3, 3),  (-6,  6,  6),  (-8,  8,  8) 

In:  (-2,  2,  2),  (-4, 4, 4),  (-5,  5,  5),  (-7,  7,  7) 

(GaP)4||(InP)4 

SI  1 6 


S.  Jungthawan  et  al. / Computational  Materials  Science  49  (2010)  SI  24— SI  18 


Table  2 

The  calculated  bandgap  deformation  coefficients  and  the  interaction  potential  VL  for 
the  alloy  configurations  in  Table  1.  The  values  in  the  parentheses  are  the  experimental 
values  of  the  coefficients  reported  in  Ref.  [41]. 


Configuration 

£(0)  (eV) 

a  (meV/GPa) 

p  (meV/GPa2) 

VL  (eV) 

InP 

Er  =  0.525 

74.1 

-1.71 

- 

El  =  1.362 

34.7 

-0.80 

GaP 

Ep  =  1.822 

89.7 

-1.44 

- 

El  =  1.616 

34.9 

-0.62 

Disordered 

Er=  1.173 

81.9  (92) 

-1.57  (-3.3) 

- 

Gao.5Irio.5P 

fn 

CO 

to 

34.8  (43) 

-0.71 

I 

E_  =  0.625 

57.3 

-1.47 

0.71 

II 

E_  =  0.708 

65.8 

-1.70 

0.60 

Ilia 

E_  =  0.762 

55.9 

-0.80 

0.56 

mb 

E_  =  0.657 

58.9 

-1.95 

0.68 

IV 

E_  =  0.662 

50.3 

-3.14 

0.67 

where  £gap(AlP),  £gap(GaP),  and  £gap(InP)  are  the  bandgaps  of  A1P, 
GaP,  and  InP,  respectively.  These  bandgap  values  from  EPM  calcula¬ 
tion  along  with  the  corresponding  experimental  values  are  given  in 
Ref.  [24]  (Table  2).  The  convergence  of  bandgap  range  is  studied  by 
calculating  the  possible  range  of  bandgaps  as  the  cell  size  increased. 
Because  the  possible  alloy  composition  of  each  cell-size  is  discrete 
in  the  composition  space,  a  third-degree  polynomial  function  of 
the  form 

A  £gap  =  a  +  bx  +  cy  +  dx2  +  ey2  +fxy  +  gx3  +  hy3  +  ixy 2 

+jx2y  (3) 

is  used  to  fit  the  calculated  bandgap  range  by  the  least  square  meth¬ 
od  to  generate  fine-mesh  data  points  for  the  contour  plots.  The  con¬ 
tour  plots  of  bandgap  range  for  different  Nmax  used  are  shown  in 
Fig.  2.  Note  that  the  bandgap  range  at  all  three  vertices  are  zero  be¬ 
cause  there  is  only  one  possible  configuration  at  those  points. 


3.  Results  and  discussion 

The  ordering  in  the  crystal  structure  induces  the  band  interac¬ 
tions  that  can  lower  the  bandgap.  Because  a  larger  cell  size  pro¬ 
vides  more  degree  of  freedom  of  ordering,  new  configurations 
with  lower  bandgaps  are  expected  as  the  cell  size  used  in  DEM  is 
increased.  The  lowering  of  the  minimum  bandgap  leads  to  the 
wider  possible  bandgap  range.  However,  the  range  of  the  bandgaps 
is  expected  to  converged  (within  a  reasonable  tolerance)  with  the 
cell  size.  The  calculated  bandgap  range  for  each  limited  cell  size  as 
a  function  of  alloy  composition  is  shown  in  Fig.  2.  We  can  see  that 
the  bandgap  range  is  quite  large  (up  to  0.8  eV)  for  the  compositions 
near  the  middle  of  the  composition  space  (slightly  shifted  toward 
the  low  Al  composition  side).  The  bandgap  range  is  reduced  as 
we  move  toward  the  vertices  and  goes  to  zero  at  the  vertices.  By 
visual  observation,  we  can  already  see  that  the  contour  plots  of 
the  bandgap  range  for  Nmax  =  5  and  above  are  quite  similar.  The 
polynomial  function  for  the  contour  plot  of  the  largest  cell  size 
(Nmax  =  8)  is  given  by 

A£gap  =  2.34x  +  4.12y  -  3.18x2  -  6.98y2  -  5.27xy  +  0.84x3 

+  2.85y3  +  2.36xy2  +  1 ,32x2y  (4) 

where  the  bandgap  reduction  (A£gap)  is  in  eV  and  x,  y  are  the  alloy 
compositions  of  AlxGayInr  _x_yP. 

To  qualitatively  determine  the  convergence  of  the  bandgap 
ranges  as  a  function  of  cell  size,  we  compute  the  difference  be¬ 
tween  the  bandgap  ranges  of  each  Nmax  with  those  of  the  largest 
cell  used  (Nmax  =  8).  This  is  done  by  calculating  the  average  root- 
mean-square  of  the  difference  (A).  The  average  is  done  over  the  en¬ 
tire  composition  space.  First,  a  200  x  200  uniform  mesh  in  the 
composition  space  is  created.  Then,  at  each  point  in  alloying  com¬ 
position  (total  of  20,301  points),  the  bandgap  range  for  each  Nmax  is 
calculated  from  the  polynomial  fitted  function.  The  root-mean- 


Fig.  2.  The  contour  plots  of  the  calculated  bandgap  ranges  of  pseudoternary  AlGalnP  as  a  function  of  alloy  compositions  (u  =  x  +  y/2  and  v  =  \f'3yf2 ).  Each  figure  shows  the 
results  of  the  calculation  with  a  different  JVmax.  Note  that  a  third-degree  polynomial  function  is  used  to  fit  to  the  calculated  data,  in  order  to  generate  the  fine-mesh  data  for  the 
plots. 
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Fig.  3.  The  root-mean-square  difference  (A)  of  bandgap  range  for  each  cell  sizes 
(Nmax).  with  respect  to  Nmax  =  8.  The  curve  represents  the  fit  to  the  data. 


£(P)  =  E(0)  +  aP  +  pP2.  The  parameters  a  and  f)  are  obtained  by  the 
least  square  fit  to  the  calculated  (first-principles)  results  and  tabu¬ 
lated  in  Table  2.  The  ordering  in  the  [111]  direction  causes  the 
conduction  band  minimum  at  L  to  fold  to  r  (called  Ecl).  The  nonlin¬ 
earity  of  £(P)  is  a  result  of  the  repulsion  between  the  conduction 
states  at  r  and  £cl  [40].  Because  the  conduction  state  is  pushed 
down,  the  bandgap  in  this  configuration  is  reduced  [25,26].  To 
study  the  effects  of  the  interactions,  the  bandgaps  of  the  ordered 
structures  are  compared  with  that  obtained  from  the  Vegard’s  rule 
(the  linear  interpolation  between  the  bandgaps  of  parent  com¬ 
pounds)  which  represents  the  case  of  completely  disordered  alloy. 
In  Table  2,  the  bandgap  of  the  disordered  alloy  is  calculated  by 
averaging  the  bandgaps  of  parent  compounds  (GaP  and  InP).  The 
interaction  between  r  and  £cl  can  be  described  by  the  first  order 
perturbation  theory.  The  perturbed  energies  of  r  and  Lc i  can  be  ex¬ 
pressed  as  [41], 


E±  =  1/2 


(EL  +  Er)±y/(EL-Er)2+4V\ 


(5) 


square  of  the  difference  (between  that  of  interested  Nmax  and  that 
of  Nmax  =  8)  in  bandgap  range  is  calculated  and  the  average  (A) 
over  all  mesh  points  are  calculated.  In  Fig.  3,  A  as  a  function  of  Nmax 
along  with  the  function  (in  the  form  of  A  =  a  +  bNc )  is  shown.  It  is 
found  that  A  started  to  converge  (with  the  tolerant  of  ~0.1  eV  ref¬ 
erenced  to  the  Nmax  =  8  calculation)  already  at  the  cell  size  contain¬ 
ing  5  cation  atoms. 

For  pseudobinary  alloys,  AIGaP  tends  to  have  a  rather  large 
bandgap  range  of  ~0.5  eV  (see  Fig.  2  for  Nmax  =  8).  AllnP  has  an 
even  larger  bandgap  range  of  ~0.6  eV.  For  GalnP,  it  is  interesting 
that  the  large  bandgap  range  of  ~0.6  eV  is  already  observed  at  a 
small  cell  Nmax  =  3.  This  indicates  that  an  ordering  that  can  signif¬ 
icantly  reduce  the  bandgap  of  GalnP  alloy  can  be  represented  by  a 
small  cell  size.  It  has  been  found  that  the  ordering  in  [111]  direc¬ 
tion  of  the  CuPt  ordered  structure  can  give  a  large  bandgap  reduc¬ 
tion  [9,24].  The  crystal  parameters  of  CuPt  ordered  structure 
consisting  of  2  cations/cell  is  shown  in  Table  1.  In  general,  we 
can  expect  more  degree  of  freedom  in  ordering  to  give  an  increased 
bandgap  range.  As  a  result,  the  smaller  bandgap  configuration 
would  occur  in  the  large  cell  size.  However,  for  this  special  case, 
the  bandgap  of  GalnP  reduction  of  CuPt  ordered  structure  provides 
the  minimum  bandgap  value  for  the  50%  alloy  (Gao.5In0.5P).  This  oc¬ 
curs  because  the  conduction  state  that  folded  onto  the  r  point 
(from  the  L  point)  interacts  strongly  with  the  conduction  band 
minimum  at  the  r  point.  Because  the  folding  that  causes  this  inter¬ 
action  already  took  place  in  the  small  cell,  increasing  cell  size  (al¬ 
low  more  flexibility  of  orderings)  does  not  make  the  interaction 
stronger  for  this  particular  case.  As  a  result,  the  bandgap  range  at 
this  composition  converges  at  a  small  cell  size.  To  further  investi¬ 
gate  the  band  interaction  for  this  particular  case,  first-principles 
calculations  are  performed. 

We  employed  first-principles  calculations  to  examine  the  elec¬ 
tronic  structures  of  several  selected  configurations  ordering  in  the 
[111]  direction.  The  crystal  parameters  and  atomic  coordinates  of 
these  configurations  are  shown  in  Table  1.  The  calculations  were 
performed  by  the  ultrasoft  pseudopotential  approach  with  the  pro¬ 
jector  augmented  wave  (PAW)  method  [34,35]  and  a  plane  wave 
basis  set  (324  eV  cutoff)  as  implemented  in  the  VASP  code  [36- 
38].  The  /"-centered  Monkhorst-Pack  k-mesh  is  used  for  Brillouin 
zone  integration. 

To  study  the  bandgap  change  as  a  function  of  pressure,  for  each 
of  alloy  configurations,  bandgap  at  the  r  point  is  calculated  at  dif¬ 
ferent  volume  (V).  To  relate  the  volume  to  the  applied  pressure,  the 
total  energy  as  a  function  of  volume  is  fitted  with  the  Rose  equa¬ 
tion  of  states  [39].  For  each  configuration,  the  pressure  dependence 
of  bandgap  (at  r)  can  be  described  by  the  relationship 


where  £+  and  £_  are  the  energies  of  perturbed  Lcl  and  r  states.  £L 
and  £r  are  the  energies  of  unperturbed  r  and  L  (obtained  by  Ve¬ 
gard’s  rule).  If  we  assumed  that  the  interaction  potential  ( VL )  does 
not  depend  on  the  pressure,  VL  can  be  calculated  by  fitting  each  sets 
of  £_,  El,  and  Er  to  Eq.  (5).  The  values  of  VL  obtained  by  the  least 
squares  fit  are  shown  in  Table  2.  The  obtained  values  of  VL  clearly 
show  that  GaP/InP  monolayer  superlattice  (structure  I  with  N  =  2) 
has  the  strongest  interaction  ( VL  =  0.71  eV).  As  a  result,  the  configu¬ 
ration  has  the  largest  bandgap  reduction  among  GalnP  alloys  [42]. 

4.  Conclusions 

The  convergence  of  the  bandgap  range  of  AlGalnP  alloys  with 
respect  to  the  cell  size  used  in  the  direct  enumeration  method 
(DEM)  is  studied.  By  directly  calculating  a  large  number  of  alloy 
configurations,  the  bandgap  range  as  a  function  of  alloy  composi¬ 
tions  are  determined  by  DEM  using  various  unit  cell  size 
(Nmax  =  2-8).  It  is  found  that  it  is  sufficient  to  use  the  cell  size  with 
5  cation  atoms  (Nmax  =  5).  The  bandgap  range  of  AlGalnP  alloy  is  re¬ 
ported  (in  the  form  of  a  polynomial  function).  It  is  found  that  the 
bandgap  range  is  large  (~0.8  eV)  for  the  compositions  near  the 
middle  of  the  composition  space  (slightly  shifted  toward  the  low 
Al  composition  side).  In  addition  to  the  convergence  study,  the 
source  of  the  bandgap  reduction  of  the  alloy  in  the  [111]  superlat¬ 
tice  configurations  are  investigated.  It  is  found  that  the  superlattice 
in  this  direction  has  strong  interactions  between  the  r  and  the 
folded  L  states  and  the  GaP/InP  monolayer  superlattice  has  the 
strong  interaction  among  studied  configurations. 

Acknowledgements 

The  work  is  supported  by  NANOTEC  (Grant  No.  NN-B-22-DI2- 
20-51-09),  AOARD/AFOSR  (Contract  No.  FA2386-09-1-4106),  and 
DOE  (Contract  No.  DE-AC36-99G010337).  The  authors  thank  L.- 
W.  Wang  for  providing  the  PESCAN  codes.  One  of  the  authors  (SJ) 
acknowledges  the  funding  from  the  Thailand  Research  Fund  (Grant 
No.  MRG5280235). 

References 

[lj  A.S.  Dissanayake,  S.X.  Huang,  H.X.  Jiang,  J.Y.  Lin,  Phys.  Rev.  B  44  (1991)  13343. 

[2]  H.X.  Jiang,  G.  Brown,  J.Y.  Lin,  J.  Appl.  Phys.  69  (1991)  6701. 

[3]  S.  Krishnamurthy,  A.  Sher,  M.  Madou,  A.B.  Chen,  J.  Appl.  Phys.  64  (1988)  1530. 

[4]  M.F.  Ling,  D.J.  Miller,  Phys.  Rev.  B  38  (1988)  6113. 

[5]  L.  Vegard,  Z.  Phys.  5  (1921)  17. 

[6]  L.  Vegard,  Z.  Kristallogr.  67  (1928)  239. 

[7]  A.  Ebina,  Y.  Sato,  T.  Takahashi,  Phys.  Rev.  Lett.  32  (1974)  1366. 

[8]  A.  Franceschetti,  A.  Zunger,  Nature  402  (1999)  60. 


SI  1 8 


S.  Jungthawan  et  al./ Computational  Materials  Science  49  (2010)  SI  14-SI  18 


[9]  K.  Kim,  P.A.  Graf,  W.B.  Jones,  J.  Comput.  Phys.  208  (2005)  735. 

[10]  L.G.  Ferreira,  S.H.  Wei,  A.  Zunger,  Int.  J.  High  Perform.  Comput.  Appl.  5  (1991) 
34. 

[11]  P.A.  Graf,  K.  Kim,  W.B.  Jones,  G.L.W.  Hart,  Appl.  Phys.  Lett.  87  (2005)  243111. 

[12]  J.  Rutherford,  Acta  Crystallogr.  A  51  (1995)  672. 

[13]  J.  Rutherford,  Acta  Crystallogr.  A  48  (1992)  500. 

[14]  L.-W.  Wang,  A.  Zunger,  J.  Chem.  Phys.  100  (1994)  2394. 

[15]  A.  Canning,  L.-W.  Wang,  A.  Williamson,  A.  Zunger,  J.  Comput.  Phys.  160  (2000) 
29. 

[16]  T.  Mattila,  L.W.  Wang,  A.  Zunger,  Phys.  Rev.  B  59  (1999)  15270. 

[17]  K.A.  Mader,  A.  Zunger,  Phys.  Rev.  B  50  (1994)  17393. 

[18]  G.L.W.  Hart,  R.W.  Forcade,  Phys.  Rev.  B  77  (2008)  224115. 

[19]  G.L.W.  Hart,  R.W.  Forcade,  Phys.  Rev.  B  80  (2009)  014120. 

[20]  K.A.  Mader,  A.  Zunger,  Phys.  Rev.  B  51  (1995)  10462. 

[21]  K.A.  Mader,  L.-W.  Wang,  A.  Zunger,  Phys.  Rev.  Lett.  74  (1995)  2555. 

[22]  A.  Franceschetti,  A.  Zunger,  Phys.  Rev.  B  52  (1995)  14664. 

[23]  J.  Kim,  L.-W.  Wang,  A.  Zunger,  Phys.  Rev.  B  56  (1997)  R15541. 

[24]  S.  Jungthawan,  S.  Limpijumnong,  R.  Collins,  K.  Kim,  P.A.  Graf  J.A.  Turner,  J. 
Appl.  Phys.  105  (2009)  123531. 

[25]  S.H.  Wei,  A.  Zunger,  Appl.  Phys.  Lett.  56  (1990)  662. 


[26]  L.C.  Su,  I.H.  Ho,  N.  Kobayashi,  G.B.  Stringfellow,  J.  Cryst.  Growth  145  (1994) 
140. 

[27]  W.  Feller,  An  Introduction  to  Probability  Theory  and  its  Applications,  vol.  1, 
Wiley,  New  York,  1968. 

[28]  A.  vandeWalle,  M.  Asta,  G.  Ceder,  CALPHAD  J.  26  (2002)  539. 

[29]  A.  van  de  Walle,  G.  Ceder,  J.  Phase  Equilib.  Diff.  23  (2002)  348. 

[30]  P.N.  Keating,  Phys.  Rev.  145  (1966)  637. 

[31]  R.M.  Martin,  Phys.  Rev.  B  1  (1970)  4005. 

[32]  J.L.  Martins,  A.  Zunger,  Phys.  Rev.  B  30  (1984)  6217. 

[33]  M.C.  Schabel,  J.L.  Martins,  Phys.  Rev.  B  43  (1991)  11873. 

[34]  P.E.  Blochl,  Phys.  Rev.  B  50  (1994)  17953. 

[35]  G.  Kresse,  D.  Joubert,  Phys.  Rev.  B  59  (1999)  1758. 

[36]  G.  Kresse,  J.  Furthmiiller,  Comput.  Mater.  Sci.  6  (1996)  15. 

[37]  G.  Kresse,  J.  Furthmiiller,  Phys.  Rev.  B  54  (1996)  11169. 

[38]  G.  Kresse,  J.  Hafner,  J.  Phys.:  Condens.  Matter  6  (1994)  8245. 

[39]  J.H.  Rose,  J.R.  Smith,  F.  Guinea,  J.  Ferrante,  Phys.  Rev.  B  29  (1984)  2963. 

[40]  T.  Kobayashi,  M.  Ohtsuji,  S.  Deol,  J.  Appl.  Phys.  74  (1993)  2752. 

[41]  J.  Dong,  Z.  Wang,  D.  Lu,  X.  Liu,  X.  Li,  D.  Sun,  Z.  Wang,  M.  Kong,  G.  Li,  Appl.  Phys. 
Lett.  68  (1996)  1711. 

[42]  S.R.  Kurtz,  J.  Appl.  Phys.  74  (1993)  4130. 


PHYSICAL  REVIEW  B  80,  193202  (2009) 


Hydrogen  doping  in  indium  oxide:  An  ah  initio  study 

Sukit  Limpijumnong,1,2,3  Pakpoom  Reunchan, 12  Anderson  Janotti,1  and  Chris  G.  Van  de  Walle1 
1  Materials  Department,  University  of  California,  Santa  Barbara,  California  93106-5050,  USA 
2School  of  Physics,  Suranaree  University  of  Technology  and  Synchrotron  Light  Research  Institute,  Nakhon  Ratchasima  30000,  Thailand 
3 Thailand  Center  of  Excellence  in  Physics  (ThEP  Center ),  Commission  on  Higher  Education,  Bangkok  10400,  Thailand 
(Received  14  August  2009;  revised  manuscript  received  13  October  2009;  published  9  November  2009) 

First-principles  density-functional  theory  is  employed  to  investigate  the  role  of  hydrogen  impurities  in 
ln203.  We  find  that  both  interstitial  hydrogen  (H;)  and  substitutional  hydrogen  (H0)  act  as  shallow  donors.  Our 
results  support  recent  experiments  by  Koida  et  al.  [Ipn.  J.  Appl.  Phys.  46,  L685  (2007)],  which  found 
hydrogen-doped  ln203  to  be  a  good  candidate  for  transparent  conducting  films. 
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Indium  oxide  (ln203)  is  one  of  the  few  transparent  mate¬ 
rials  that  can  be  doped  to  a  very  high  carrier  concentration 
without  significantly  degrading  their  transparencies.  The 
most  widely  used  transparent  conductor  is  tin-doped  ln203, 
known  as  indium  tin  oxide  (ITO),  which  is  used  as  a  trans¬ 
parent  electrode  for  optoelectronic  devices  such  as  solar  cells 
and  light  emitting  diodes.  Recently,  Koida  et  al.1  showed  that 
hydrogen-doped  ln203  can  have  mobility  exceeding 
100  cm2/V  s  at  carrier  densities  of  ~  1020  cm-3  and  exhibits 
good  transparency  even  at  near-infrared  (NIR)  wavelengths. 
Hydrogen-doped  ln203  is,  therefore,  a  promising  material  to 
replace  or  improve  upon  ITO  for  applications  in  optoelec¬ 
tronic  devices  such  as  solar  cells  or  photodetectors,  in  par¬ 
ticular,  when  sensitivity  in  the  NIR  region  is  required. 

In  this  Brief  Report,  we  provide  a  detailed  first-principles 
study  of  H  in  ln203.  In  most  semiconductors  H,  is  amphot¬ 
eric,  acting  as  a  donor  in  p-type  samples  and  as  an  acceptor 
in  u-type  samples.  If  H,  is  amphoteric,  it  can  never  be  the 
cause  of  conductivity  since  it  self-compensates.  There  are 
exceptions  though.  For  example,  H,  was  proposed,  based  on 
first-principles  calculations,2  to  act  exclusively  as  a  donor  in 
ZnO,  a  prediction  which  was  soon  confirmed  by 
experiments.^'  Hydrogen  is  thus  one  of  the  causes  of  uninten¬ 
tional  ?;-type  conductivity  in  ZnO.  In  the  present  work  we 
find  that  H,  is  also  a  shallow  donor  in  ln203,  and  can,  there¬ 
fore,  also  be  the  cause  of  n-type  conductivity.  In  addition,  H 
atoms  can  also  occupy  substitutional  oxygen  sites  (H0)  and 
also  act  exclusively  as  donors  in  this  configuration.  Substitu¬ 
tional  hydrogen  was  recently  described  in  terms  of  a  multi¬ 
center  bond  and  found  to  be  stable  in  a  number  of  oxides  and 
some  nitride  materials.4"6  In  order  to  explore  the  stability  of 
H0,  we  also  performed  calculations  for  oxygen  vacancies 
(V0)  in  ln203.  Our  findings  confirm  that  hydrogen  acts  as  the 
source  of  u-type  conductivity  in  ln203. 

In203  has  a  fairly  complicated  crystal  structure  called  bix- 
byite  [space  group  T7h(Ia3 )]  with  80  atoms  in  a  conventional 
cell  and  40  atoms  in  the  primitive  unit  cell.  Its  large  primitive 
unit  cell,  compared  with  conventional  semiconductors  (two 
atoms  per  cell  for  Si  or  GaAs),  rendered  first-principles  cal¬ 
culations  of  this  material  prohibitively  demanding  in  the 
past.  Only  in  the  last  decade  have  first-principles  calculations 
for  ln203  been  attempted,  starting  with  a  linear  muffin  tin 
orbital  method  with  atomic  sphere  approximations  with  a 
small  basis  set  for  an  unrelaxed  experimental  structure  by 


Odaka  et  al.,1  followed  by  a  more  advanced  full-potential 
linear  muffin-tin-orbital  calculation  that  included  volume  re¬ 
laxation  by  Mryasov  and  Freeman.8  Based  on  their  investi¬ 
gation  of  the  band  structure,  Mryasov  and  Freeman  attributed 
the  coexistence  of  conductivity  and  transparency  to  the 
unique  characteristics  of  the  conduction-band  separation  in 
ln203.  More  recently,  Medvedeva9  used  the  full-potential  lin¬ 
earized  augmented  plane-wave  method  and  allowed  full 
structural  relaxations  to  calculate  the  band  structures.  To 
date,  only  a  few  impurities  in  ln203  have  been  studied  by 
first-principles  calculations.8"11 

We  perform  first-principles  calculations  using  density- 
functional  theory  (DFT)  within  the  local-density  approxima¬ 
tion  (LDA),  as  implemented  in  the  Vienna  ab  initio  simula¬ 
tion  package.12  Electron-ion  interactions  are  treated  using 
projector  augmented  wave  potentials.13  The  electron  wave 
functions  are  described  using  a  plane-wave  basis  set  with  an 
energy  cutoff  of  400  eV.  Since  In  4 d  electrons  are  not  fully 
localized  near  the  core  and  can  play  a  role  in  the  crystal 
bonding,  they  are  treated  as  valence  electrons  and  allowed  to 
relax  during  the  self-consistent  calculations.  A  shifted  2X2 
X  2  /.'-point  sampling  based  on  the  Monkhorst-Pack  scheme14 
is  employed  for  the  Brillouin-zone  integration.  The  calcu¬ 
lated  lattice  constant  of  10.073  A  is  in  good  agreement  with 
the  experimental  value  of  10.117  A.15  The  calculated  band 
gap  (1.17  eV)  is  much  smaller  than  the  experimental  value 
due  to  the  well-known  DFT  problem.  In  the  past,  the  experi¬ 
mental  band  gap  of  ln203  was  widely  believed  to  be  3.6 
eV. 16  Only  recently  has  the  actual  value  of  the  band  gap  been 
established,  at  2.67  eV  for  molecular-beam  epitaxy  (MBE)- 
grown  samples17  and  a  slightly  higher  value  at  2.89  eV  for  rf 
magnetron  sputtered  samples.18  We  can  overcome  the  band- 
gap  underestimation  inherent  in  LDA  by  performing  LDA 
+  U  calculations,  allowing  us  to  investigate  the  changes  in 
the  transition  levels  and  formation  energies  as  the  band  gap 
is  increased.  This  approach  has  been  successfully  applied  in 
the  study  of  H  in  other  oxides  and  nitrides.4"6,19  To  calculate 
H  in  various  configurations,  we  follow  a  supercell  approach 
with  a  conventional  unit  cell  of  80  atoms  as  the  supercell  and 
the  lattice  constant  fixed  at  the  calculated  bulk  value 
(10.073  A).  All  atoms  are  allowed  to  relax  during  the  calcu¬ 
lations.  The  formation  energy  of  Hf  {q=—,  0,  or  +)  is  defined 
as20,21 
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E'X Hf)  =  £tot(Hf)  -  £tot(bulk)  -Mh  +  qEF,  (1) 

where  £tot(H?)  is  the  total  energy  of  the  supercell  containing 
H,  in  charge  state  q,  and  £’tol(bulk)  is  the  total  energy  of  the 
perfect  crystal  in  the  same  supercell.  For  a  charged  H,  (q 
#0),  electrons  are  exchanged  with  a  reservoir  with  energy 
EF,  the  Fermi  energy  of  the  system.  /rH  is  the  hydrogen 
chemical  potential,  i.e.,  the  energy  of  the  reservoir  with 
which  H  atoms  are  exchanged.  We  fixed  the  value  of  fiu  to 
half  of  the  calculated  energy  of  a  free  H2  molecule.  The 
formation  energies  of  Vq  and  H0  are  defined  as 

Ef(Vq0)  =  EUK)  -  £to,(bulk)  +  /r0  +  qEF,  (2) 

and 

Ef (Ho)  =  £tot(HS)  -  £,ot(bulk)  +  /r0  -  /rH  +  qEF,  (3) 

where  E^V^)  is  the  total  energy  of  the  supercell  with  one  O 
removed,  and  £tot(H^)  is  the  total  energy  of  the  supercell 
with  one  H  substituted  for  O.  For  the  purposes  of  presenting 
our  results  we  assumed  In-rich  growth  conditions,  i.e.,  /rIn  is 
fixed  at  the  energy  per  atom  of  In  metal,  and  fi0  is  fixed  at 
3.66  eV  below  the  energy  of  half  an  02  molecule,  based  on 
the  calculated  heat  of  formation  of  ln203  (10.99  eV/f.u.). 
Values  for  other  growth  conditions  can  easily  be  obtained  by 
referring  back  to  Eqs.  (2)  and  (3). 

By  symmetry,  ln203  in  the  bixbyite  structure  has  two  in¬ 
equivalent  In  sites  (8 b  and  24 d  in  Wyckoff  notation)  and  all 
oxygen  sites  are  equivalent.1'’  A  quarter  of  all  In  atoms  (la¬ 
beled  Ini)  occupy  the  8 b  positions  while  the  remaining  three 
quarters  (labeled  In2)  occupy  the  2Ad  positions.  Both  types 
of  In  atoms  have  a  sixfold  coordination  with  their  O  neigh¬ 
bors.  The  calculated  local  structures  of  Ini  and  In2  atoms  are 
shown  in  Figs.  1(a)  and  1(b),  respectively.  The  local  structure 
of  Ini  atom  is  quite  symmetrical.  All  six  In-O  bonds  are 
equivalent  with  a  bond  distance  of  0.2%  smaller  than  the 
average  value  of  all  In-O  bonds.  It  has  a  shear  distorted 
octahedral  coordination  as  shown  in  Fig.  1(a).  If  all  the  bond 
angles  were  90°,  the  structure  would  be  a  perfect  octahedron. 
The  In2  configuration  is  less  symmetrical,  as  shown  in  Fig. 
1(b).  For  In2,  the  six  In-O  bonds  can  be  divided  into  three 
groups  of  two  equivalent  In-O  bonds.  The  major  difference 
between  Ini  and  In2  is  that  Ini  and  its  O  neighbors  form 
three  sets  of  linear  O-In-O  bonds  whereas  no  linear  O-In-O 
bonds  occur  for  In2.  For  the  local  structure  of  oxygen,  each 
O  atom  forms  distorted  tetrahedral  bonds  with  its  neighbor¬ 
ing  In  atoms  (one  Ini  and  three  In2)  as  shown  in  Fig.  1(c). 
These  four  In-O  bonds  are  inequivalent  and  have  bond  dis¬ 
tances  that  differ  from  the  average  bond  distance  ( 

=  2.169  A)  by  +2.0%,  +0.5%,  -0.2%,  and  -2.3%. 

We  have  calculated  H,  in  three  charge  states:  +,  neutral, 
and  -.  Similar  to  other  oxides2'22  and  nitrides,20,23  H3-  prefers 
sites  near  the  anions.  In  other  oxides  and  nitrides,  H3  exhibits 
local  minima  at  the  bond  center  (BC)  and  anion  antibonding 
(AB0)  sites.  For  ln203  there  are  four  inequivalent  BC  and 
four  inequivalent  AB0  sites,  as  illustrated  in  Fig.  1(d).  We 
find  that  H3  is  not  stable  at  the  BC  sites.  It  is  spontaneously 
pushed  away  from  the  site  and  moves  significantly  off-axis 
(similar  to  the  OAy  site  in  Ref.  24).  Indeed,  two  of  the  BC 
sites  lead  to  configurations  that  are  essentially  the  same  as 


FIG.  1 .  (Color  online)  Schematic  illustration  of  the  coordination 
around  In  and  O  atoms  in  ln203.  The  percentage  numbers  indicate 
the  differences  in  bond  length  compared  to  the  average  of  all  In-O 
bonds  (calculated  value  =  2.169  A).  There  are  two  In  inequivalent 
sites,  both  surrounded  by  six  O  atoms.  All  O  sites  are  equivalent 
and  surrounded  by  four  In  atoms,  (a)  The  8 b  In  site,  (b)  the  24 d  In 
site,  (c)  the  O  site,  and  (d)  the  possible  sites  for  interstitial  H. 

the  AB03  and  AB04  sites.  The  calculated  energies  for  H~  at 
the  four  inequivalent  AB0  sites  are  listed  in  Table  I.  It  has 
been  previously  shown  for  the  cases  of  ZnO  (Ref.  2)  and 
Ill-Nitrides20,23  that  when  H  occupies  the  AB  site,  the  corre¬ 
sponding  host-atom  bond  is  significantly  extended.  We  ob¬ 
serve  the  same  phenomenon  in  ln203.  This  may  explain  why 

TABLE  I.  Calculated  formation  energies  E f  for  H  defects  and  O 
vacancies  in  ln203.  For  some  defects  with  an  occupied  gap  state,  a 
corrected  value  of  E?  based  on  LDA+  U  and  extrapolation  (Ref.  19) 
(see  text)  is  presented  in  square  brackets.  For  charged  defects  we 
take  the  Fermi  energy  at  the  VBM.  Equilibrium  with  bulk  In  metal 
(In-rich  conditions)  and  with  H2  molecules  at  T=  0  is  assumed. 
Bond  lengths  for  O-H  bonds  in  the  antibonding  configurations  of 
H3"  are  also  listed. 


Defect 

Location 

Ef 

(eV) 

4o-h 

(A) 

nr 

AB0i 

-2.09 

1.026 

nr 

abG2 

-1.66 

1.015 

nr 

ab03 

-1.34 

0.995 

»r 

ab04 

-1.35 

0.996 

H° 

16c 

2.43 

nr 

16c 

3.41  [4.27] 

nr 

8  a 

4.33 

hs 

O  site 

-1.40 

Vo 

O  site 

-2.68 

k 

O  site 

-0.84[0.01] 

K 

O  site 

0.96  [2.30] 
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FIG.  2.  (Color  online)  Schematic  illustration  of  H  defects  in  ln203.  Hi1"  at  the  AB01  site  in  (a),  H(  at  the  16c  site  in  (b),  and  Hq  in  (c). 


the  AB01  site  is  the  preferred  site  for  Hi1",  since  its  corre¬ 
sponding  In-O  bond  is  initially  the  longest.  Ht  at  the  AB01 
site  is  lower  in  energy  than  the  other  (metastable)  sites  by  at 
least  0.4  eV,  rendering  the  occupation  of  those  other  sites 
negligible  at  reasonable  temperatures.  The  fully  relaxed  ge¬ 
ometry  of  H7"  at  the  AB01  site  is  shown  in  Fig.  2(a).  The  O-H 
bond  distances  (r/0_H),  as  shown  in  Table  I,  are  slightly 
longer  than  the  O-H  bond  in  HiO.  A  similar  trend  was  ob¬ 
served  for  the  case  of  H,  in  ZnO.27 

In  the  negative  charge  state,  H,  prefers  interstitial  sites 
closer  to  the  cations.20,23  In  the  bixbyite  crystal,  two  such 
interstitial  sites  exist,  which  we  will  refer  to  as  8a  and  16c 
following  their  Wyckoff  notations. 15  We  find  that  Hf  and  H° 
are  stable  at  the  16c  site  with  the  single-electron  level  at 
~1  eV  above  the  valence-band  maximum  (VBM);  the  cor¬ 
responding  defect  state  is  spatially  localized  on  the  H  atom. 
The  relaxed  geometry  of  H~  at  the  16c  site  is  shown  in  Fig. 
2(b).  At  the  8a  site  Hf  is  metastable  with  a  formation  energy 
that  is  higher  by  0.9  eV  than  at  the  16c  site. 

The  formation  energies  for  H“  and  H°  are  shown  in  Table 
I,  and  a  plot  of  formation  energy  as  a  function  of  Fermi-level 
position  for  all  three  charge  states  is  shown  in  Fig.  3.  Based 
on  the  LDA  results,  H,  is  a  negative-f/  defect  (the  neutral 


FIG.  3.  (Color  online)  Calculated  formation  energies  of  hydro¬ 
gen  defects  and  oxygen  vacancies  in  ln203  as  a  function  of  Fermi 
energy.  In-rich  conditions  are  assumed. 


charge  state  is  never  stable).  H,  is  stable  in  the  +  charge  state 
for  the  range  of  Fermi  energies  from  the  VBM  up  to  EF 
~2.75  eV  above  the  VBM.  Since  the  band  gap  of  MBE- 
grown  ln203  has  been  reported  to  be  2.67  eV,17  we  can  con¬ 
clude  that  hydrogen  behaves  exclusively  as  a  donor  for 
Fermi  levels  within  this  gap. 

The  bare  LDA  results  show  that  the  +/-  transition  level 
(the  Fermi-level  value  for  which  Hi1"  and  Hi"  have  equal  for¬ 
mation  energies)  occurs  at  2.75  eV  above  the  VBM.  Since 
this  value  is  only  slightly  above  the  reported  gap  of  2.67 
eV,17  and  since  there  may  still  be  some  uncertainty  in  this 
band-gap  value  (due  to  the  band-gap  renormalization  and  the 
potential  residual  strain  in  the  thin  films),  it  is  worthwhile  to 
obtain  a  more  reliable  value  for  the  +/-  transition  level  by 
going  beyond  the  bare  LDA  results.  LDA  underestimates  the 
band  gap  and,  therefore,  also  the  energy  of  Hf,  which  has  an 
occupied  state  in  the  band  gap.  Correction  of  the  band  gap 
may  shift  this  state  upward  and  hence  increase  the  energy  of 
Hf.  To  estimate  this  correction,  we  performed  LDA+f/  cal¬ 
culations  as  described  in  Ref.  19  with  U=4  eV  for  the  In  4 d 
states.  The  LDA+t/  calculation  for  bulk  ln203  produces  a 
slightly  reduced  lattice  constant  (by  1.4%)  and  raises  the 
band  gap  by  0.47  eV  compared  to  LDA.  To  obtain  a  com¬ 
plete  band-gap  correction,  further  extrapolation  is  needed.19 
Based  on  the  calculations  of  bulk  band  gap,  we  can  see  that 
the  correction  has  to  be  scaled  up  by  (2.67- 1.17) /0.47 
=  3.19  times  to  get  the  full  band-gap  correction.  The  LDA 
+  U  calculations  of  Hf  raise  the  +/-  level  by  0.14  eV;  there¬ 
fore,  the  corrected  transition  level  would  occur  at  2.75 
+  0.14X3.19=3.20  eV,  which  is  well  above  the  conduction- 
band  minimum  (CBM).  The  corrected  formation  energy  of 
H  7  is  shown  in  Lig.  3. 

In  addition  to  H„  we  have  studied  H  substitution  on  the  O 
site,  finding  that  H0  is  also  exclusively  a  donor.  It  always 
exists  in  +  charge  state  without  any  defect  level  inside  the 
band  gap.  The  fully  relaxed  geometry  of  Hq  is  shown  in  Lig. 
2(c),  and  the  calculated  formation  energy  of  H0  is  included 
in  Table  I  and  plotted  in  Lig.  3.  The  formation  energy  of  H0 
is  slightly  higher  than  that  of  H„  but  low  enough  to  lead  to  a 
significant  incorporation  of  substitutional  hydrogen.  We  also 
performed  calculations  for  oxygen  vacancies  in  order  to  in¬ 
vestigate  the  stability  of  H0  with  respect  to  the  dissociation 
into  H,  and  VQ.  We  find  that  VQ  is  a  deep  donor  in  ln203.  Our 
calculations  show  that  V0  is  stable  in  two  charge  states:  2+ 
and  neutral;  the  +  charge  state  is  never  stable.  We  again  used 
the  LDA  +  U  extrapolation  scheme;  formation  energies  based 
on  LDA  and  LDA+  U  are  listed  in  Table  I,  and  the  corrected 


193202-3 


BRIEF  REPORTS 


PHYSICAL  REVIEW  B  80,  193202  (2009) 


formation  energies  are  shown  in  Fig.  3.  Our  calculated  2 
+  /  0  level,  referenced  to  the  VBM,  is  in  good  agreement  with 
that  calculated  by  Lany  and  Zunger.26  (Note  that  in  Lany  and 
Zunger’s  work,  the  transition  level  for  V0  appears  to  be 
much  deeper  in  the  gap  because  they  used  a  larger  ln203 
bandgap  of  3.50  eV.)  The  binding  energy  between  V()  and  H, 
to  form  H0  can  be  calculated  from  the  formation  energies  of 
the  three  defects,  shown  in  Fig.  3.  Under  u-type  conditions 
the  binding  energy  is  1.82  eV,  a  large  value  which  is  indica¬ 
tive  of  the  stability  of  the  substitutional  hydrogen. 

Recently,  Koida  et  al. 1  prepared  and  characterized  rf  mag¬ 
netron  sputtered  ln203  samples  doped  with  various  concen¬ 
trations  of  hydrogen.  They  reported  an  optimized  carrier  den¬ 
sity  of  3  X  1020  cm-3  and  mobility  as  high  as  130  cm2/V  s 
for  a  sample  doped  with  ~3  at.  %  of  H.  King  et  al.21  per¬ 
formed  muon  spin  rotation  measurements  to  investigate  the 


electrical  behavior  of  hydrogen  in  ln203.  Their  results  indi¬ 
cate  that  hydrogen  is  a  shallow  donor,  with  the  estimated 
+/-  transition  level  at  0.6  eV  above  the  conduction-band 
minimum.27  These  experimental  results  are  consistent  with 
our  calculations,  which  indicate  that  both  interstitial  and  sub¬ 
stitutional  hydrogen  serve  as  donors  in  ln203. 
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